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Samenvatting 



De komst van ultrasnelle computers heeft een enorme vooruitgang betekend voor 
de simulatie van kwantumveeldeeltjesproblemen. Ze zijn een essentieel onderdeel 
geworden van de studie van complexe systemen gaande van biologische toepassin- 
gen tot de materiaalwetenschappen. Deze systemen vereisen allemaal de numerieke 
oplossing van de elektronische Schrodinger-vergelijking zodanig dat de totale ener- 
gie, de elektronendichtheid of andere eigenschappen van het systeem voor handen 
zijn. Wanneer de Schrodinger-vergelijking in zijn exacte vorm wordt opgelost, botst 
men al snel op de exponentiele toename van de hoeveelheid informatie die vereist 
is. Dit maakt dat een exacte oplossing onmogelijk is en er benaderingen dienen 
gemaakt te worden. 

Het doel van de kwantumchemie is de bepaling van de elektronische structuur tot op 
chemische nauwkeurigheid. Dit betekent: vergelijkbaar met experimentele waarden 
en met voorspellende kracht. Elke ontwikkelde methode moet onderworpen worden 
aan zorgvuldige screening van deze eis door middel van testsystemen. 
Hartree-Fock is een theorie die soms voldoende betrouwbare resultaten oplevert. 
Maar toch is het met deze techniek onmogelijk om correcte voorspellingen te be- 
komen voor chemische processen waarin elektroncorrelatie belangrijk wordt. De 
zoektocht naar een techniek die zowel betrouwbaar is als computationeel goed 
schaalt, gaat dus onvermoeibaar verder. 

Het ontwikkelen van nieuwe kwantummechanische methoden gebeurt deels met het 
oog op het gebruik in industrieel onderzoek door niet-specialisten. Daardoor is het 
belangrijk om ervoor te zorgen dat de methode zonder veel parameters of instel- 
lingen kan toegepast worden. Ab initio methoden hebben als voordeel dat er geen 
empirische data vereist is (en dus ook geen fundamenteel begrip) om berekingen 
te doen. Deze methoden gebruiken enkel fundamentele waarden verbonden aan 
het systeem. Dit maakt ze ideale kandidaten voor dagelijks gebruik in kwantum- 
chemische berekeningen. Het is echter niet eenvoudig een kwantummechanische 
techniek te vinden die over een wijd spectrum van systemen kan toegepast worden. 
Greense functies zijn de aangewezen grootheden voor de directe bepaling van de 
elektronische structuur van atomen en molecules. In plaats van de volledige elek- 
tronische golffunctie te berekenen, volstaat het om de eigenschappen van enkele 
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elektronen in hun omgeving te beschouwen. De overeenkomstige fysische groothe- 
den zijn de Greense functies of propagatoren. Er zijn verscheidene voordelen aan 
het gebruik van Greense functies bij de beschrijving van de fysica van veeldeel- 
tjessystemen. Ten eerste vormen ze een directe berekeningsmethode voor fysische 
grootheden. Eens de Greense functie gekend is, is het mogelijk de verwachtings- 
waarde van elke eendeeltjesoperator en de grondtoestandsenergie te bepalen zon- 
der verdere berekeningen. Ook de ionisatie-energieen en spectrale sterkte is binnen 
handbereik. Ten tweede laten ze een eenvoudige fysische interpretatie toe. De 
kwantummechanische Greense functie bepaalt de evolutie van toestanden in de 
tijd en ruimte, net als de klassieke tegenhanger. De diagrammatische representa- 
tie in Feynman-diagrammen biedt een gemakkelijke en universele mogelijkheid om 
theorieen te vergelijken en te verduidelijken. Op deze manier kan men ook nieuwe 
methoden ontwikkelen. 

In dit werk introduceren we een Greense-functietechniek voor de berekening van 
grondtoestandsenergieen en ionisiatie-energieen in kwantumveeldeeltjessystemen. 
De elementaire bouwstenen van deze theorie zijn de RPA (random phase approxi- 
mate) excitaties. Deze worden gebruikt om een benadering voor de zelfenergie op 
te bouwen via een stelsel Faddeev-vergelijkingen. De resulterende FRPA (Faddeev 
RPA) is een potentieel wijd toepasbare Greense functie techniek. Het gebruik van 
de RPA heeft als gevolg dat de techniek bruikbaar is voor uitgestrekte systemen 
zoals het uniform elektronengas en nucleaire materie, maar ook voor eindige syste- 
men zoals atomen en moleculen. Wanneer de FRPA in zijn volledig zelfconsistente 
vorm wordt uitgevoerd, voldoet deze aan de Baym-Kadanoff voorwaarden en zal 
er bijgevolg aan belangrijke behoudswetten voldaan zijn. 

Eerst worden de theoretische concepten ge'i'ntroduceerd die nodig zijn om een be- 
nadering in te voeren voor de zelfenergie. De Dyson-vergelijking wordt gefor- 
muleerd in functie van de zespunts vertexfunctie. Dit maakt de koppeling tus- 
sen eendeeltjestoestanden en 2plh (two-particle-one-hole) en 2hlp (two-hole-one- 
particle) toestanden natuurlijk. Daarna worden de polarisatiepropagator en de 
tweedeeltjes Greense functie ingevoerd. De toepassing van de RPA voor deze ob- 
jecten vormt een benadering voor de ph en pp excitaties. 

In Hoofdstuk 3 wordt de FRPA techniek afgeleid. Het resultaat is een methode 
die pp (particle-particle) en ph (particle-hole) interacties van het RPA niveau ver- 
mengt op gelijke voet. De procedure die daarvoor zorgt, is ontwikkeld in een 
ruimte die naast de eendeeltjestoestanden ook de 2plh en 2hlp toestanden om- 
vat. De Faddeev-vergelijkingen zorgen ervoor dat de pp en ph excitaties op een 
correcte manier gekoppeld worden tot 2plh en 2hlp toestanden. De eigenwaarden 
en transitie-amplitudes uit het RPA probleem worden gebruikt om de Faddeev- 
matrices op te bouwen. Het is absoluut noodzakelijk dat de Faddeev-techniek 
gebruikt wordt om excitaties te kunnen koppelen van het RPA type. Het is op 
dit vlak dat de gebruikte techniek verschilt van de gangbare ADC(3) (algebraic 
diagrammatic construction tot op derde orde); een techniek die interacties bevat 
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van het TDA (Tamm-Dancoff approximation) type. Niettegenstaande zijn ze beide 
toch erg gelijkend. De ADC(3)-vergelijkingen kunnen afgeleid worden uit de FRPA- 
vergelijkingen door het weglaten van de achterwaards propagerende amplitudes van 
het RPA-probleem. We noemen ADC(3) vanaf nu dan ook FT DA (Faddeev TDA). 
Deze equivalentie was al duidelijk uit diagrammatische overwegingen, maar wij 
hebben dit nu ook op het vlak van de matrixvergelijkingen aangetoond. Doordat 
de FRPA de TDA transitie-amplitudes en eigenwaarden vervangt door hun RPA te- 
genhangers, hoort deze beter de correlaties te beschrijven in uitgestrekte systemen. 
Op het eerste gezicht lijkt het alsof de FRPA de matrixdimensies met een factor 
drie doet toenemen. Wanneer men echter de aard van de oplossingen bekijkt, dan 
kan men besluiten dat twee derden van deze oplossingen onfysisch zijn en kunnen 
weggelaten worden. De fysische en onfysische oplossingen kunnen gemakkelijk on- 
derscheiden worden. Een matrixprojectie herleid de matrixdimensies tot die van 
een 2plh en 2hlp diagonalisatie. 

Hoofdstuk 4 bevat de resultaten die berekend zijn met de FRPA. Als eerste toepas- 
sing is de elektronische structuur van een reeks atomen berekend. Deze berekenin- 
gen tonen aan dat zelfs voor erg kleine systemen, waar de RPA het Pauli-principe 
hard schendt, de berekende waarden van goede kwaliteit zijn. Algemeen kunnen 
we zeggen dat FRPA en FTDA zeer gelijklopende resultaten geven voor de lichtste 
atomen. Naarmate het atoomgetal toeneemt, worden de additionele diagrammen 
in de RPA van groter belang en vormen de FRPA-waarden een verbetering ten 
opzichte van FTDA. Beide methoden volgen de trend van CCSD (coupled cluster 
met singles en doubles) in de limiet van oneindig grote basis set en de afwijking ten 
opzichte van het experiment ligt binnen de afgeschatte extrapolatiefout. De bijna 
ontaarding in het Be atoom vormt een groter probleem. De fouten, die bij elke 
methode voorkomen, zijn waarschijnlijk vooral te wijten aan tekortkomingen van 
de basis set. Voor de ionisatie-energieen gelden dezelfde conclusies: de kleinere 
twee-elektron atomen vertonen geen verschillen tussen FRPA en FTDA, terwijl voor 
alle andere atomen de FRPA het resultaat in de richting van het experiment duwt. 
Ook hier vormt het Be atoom een problematisch geval. 

Na deze berekeningen voor atomen is het logisch om over te gaan tot het bepalen 
van de eigenschappen van molecules. Bij molecules vormt de internucleaire afstand 
een extra parameter. Rond de evenwichtsgeometrie vinden we dat de resultaten 
gelijklopend zijn met die uit de atomaire berekeningen. Voor kleinere molecules 
vertonen FRPA en FTDA weinig verschillen, terwijl voor grotere molecules de ener- 
gieen verder uit elkaar liggen. Het aflopen van verschillende internucleaire afstan- 
den toont in elk geval aan dat de zelfconsistentie voor het Hartree-Fock diagram 
noodzakelijk is. Zonder deze iteratieve procedure wordt er voor sommige molecules 
zelfs geen minimum gevonden in de potentiele-energiecurve. De ionisatie-energieen 
zijn ongeveer even goed in FRPA als in FTDA. Er is net als bij atomen ook een 
problematisch geval bij het bepalen van de ionisatie-energieen: de derde ionisatie- 
energie van de N2 molecule wijkt sterk af van het experiment. In dit hoofdstuk 
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werd ook een poging ondernomen om de basis set limiet te bereiken, maar dit ligt 
nog buiten de computationele mogelijkheden. 

Moleculaire berekeningen in de dissociatielimiet leggen het zwakke punt van de 
FRPA bloot. Doordat de FRPA voor de opbouw van de matrices zeer sterk afhangt 
van de RPA, kan de FRPA ook enkel toegepast worden daar waar de RPA een zinvol 
resultaat produceert. Het is bekend dat de ph RPA een triplet instabiliteit vertoont 
wanneer de internucleaire afstand van de diatomische moleculen vergroot wordt. Bij 
de kleinere molecules vinden we niet dat de kwaliteit van de resultaten verslechtert 
in de buurt van de instabiliteit. Voor de grotere moleculen is dit echter niet geldig, 
zoals de berekeningen voor de HF molecule aantonen. Bij molecules met een hogere 
correlatie-energie wijkt het FRPA-resultaat sterk af van de verwachte resultaten in 
de buurt van de instabiliteit. Dit betekent dat de FRPA vooral een methode is die 
rond de evenwichtsgeometrie geldig is. 

De Hubbard-molecule is een testsysteem dat hetzelfde gedrag vertoont als een 
molecule in de dissociatielimiet. Dit model laat toe om de problematiek van disso- 
ciate te bekijken in de eenvoudigst mogelijke vorm. De exacte resultaten voor dit 
model vormen een leidraad voor het aanbrengen van verbeteringen in de huidige 
theorie. Door fragmentatie toe te laten in de eendeeltjespropagator is het mogelijk 
om voorbij de instabiliteit van de RPA toch berekeningen te doen met de FRPA. 
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Chapter 1 

Introduction 



The development of modern computers over the last 50 years has enabled the 
simulation of the quantum many-body problem[DN05, FW71, RS80]. They have 
become an essential tool for studying complex systems such as those arising in 
biology and material science. The underlying technology is the computational 
solution of the electronic Schrodinger equation, i.e. to calculate the electronic 
energy, electron density and other properties. In its exact form, this Schrodinger 
equation becomes exponentially harder to solve with increasing number of electrons. 
Hence, a brute force solution is out of the question and an approximation given 
the positions of the nuclei and the number of electrons has to be made. 
A large number of chemical properties can be calculated, among which the elec- 
tronic structure and spectroscopic constants of molecules. The ambition of the 
quantum chemist has always been to predict these properties with chemical ac- 
curacy, meaning that the results obtained are able to compete with experimental 
data. Every method that is developed needs to be submitted to careful testing and 
benchmarking on test systems to prove its value. 

Hartree-Fock theory produces results that are reasonably accurate, but is incapable 
of yielding correct results in chemical events where correlation is important. Thus, 
finding a robust treatment of electron correlation that exhibits a tractable scaling in 
computational effort with the size of the system, is the key problem for present-day 
model developers. 

The aim of a newly developed method is always the widespread use by non- 
specialists in research and industry. With this in mind, the black-box character of 
a method greatly determines its applicability in the wider field. Ab initio methods 
constitute a class of simulation techniques that do not rely on empirical input (and 
thus deeper understanding) from the user. These methods use only fundamental 
knowledge of the system of interest. This means that they are the perfect can- 
didates for every-day use in quantum chemical calculations. Nevertheless, it is not 
easy to find an ab initio method that can be used for any system from the smallest 
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toy models to extended systems such as solids. 

1.1 Motivation 

Green's function methods[DN05, FW71, LO05] represent a very useful set of 
tools for direct calculation of electronic structure properties of molecules[Danll, 
vNSC84, TF74, Kut09]. They rely on the idea that for the description of various 
properties of the molecule, it is not necessary to know the exact information of 
every electron. It suffices to describe one or two particles within their environ- 
ment. The corresponding physical quantities are the one- and two-particle Green's 
functions. 

The advantages of using Green's functions to describe physical phenomena in many- 
body systems are numerous. First of all, the most important physical quantities 
can be calculated directly. When the Green's function is known, deriving these 
variables requires no further calculations. Secondly, they have a simple physical 
interpretation. Just like the classical analogon, the quantum mechanical Green's 
function can be seen as the object that propagates the particles of interest in 
time and space. Another advantage is the easy and systematic way of expressing 
Green's functions in terms of Feynman diagrams. These graphical representations 
can be used as a tool for the interpretation of certain approximations and are a 
powerful communication tool about different theories. They also present a succesful 
construction path for new approximation schemes. 

The Green's function formalism provides a good alternative to widely spread many- 
body techniques such as quantum Monte Carlo (QMC), configuration interaction 
(CI) and coupled cluster (CC) theory. 

The remainder of this introduction presents a short overview of the historical evol- 
ution of methods as well as the relevant literature that has led to the development 
of the approach considered in the present work. 

1.1.1 Extended systems 
1.1.1a Electron gas 

A wide array of test systems has been developed to study the behavior of the Cou- 
lomb interaction in atoms, molecules and solids. The uniform electron gas consists 
of a uniform negative charge distribution that is governed by a Coulombic force 
against a backdrop of positive charge. There has been much interest pioneered 
by Pines[Pin71] in describing the response of this system to an external perturb- 
ation. The object that holds information about this response is the polarization 
propagator. Apart from a continuum spectrum of particle-hole (ph) poles, this 
polarization propagator exhibits a discrete pole in the low momentum region cor- 
responding to the plasmon excitation. When describing the polarization propagator 
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in the Tamm-Dancoff approximation, the pole energy diverges due to the beha- 
vior of the Coulombic interaction. This is illustrated in Figure 1.1. The lowest level 
of theory that describes the plasmon pole correctly is the random phase approx- 
imation (RPA)[BP51, PB52, BP53]. The RPA polarization screens the Coulomb 
interaction in such a way that the classical plasmon energy 



is correctly reproduced. Here p is the uniform density of the electron gas. This 
means that for a method to be applicable in this low momentum limit, correlations 
at least on the RPA level or beyond need to be included. 

The GW approximation (GWA)[Hed65, AG98, ORR02] applied in the electron gas 
theory is one of the methods that benefits from the succesful description of the 
screening of the interaction in RPA. This method can be seen as a version of 
Hartree-Fock where the bare interaction is replaced by a dynamically screened 
interaction. The success for the electron gas raises the question whether this 
method can also be applied to smaller systems. These systems are characterized by 
a stronger electronic relaxation effect and weaker screening than extended systems. 
Research has been conducted in continuous[RJT10] and discrete[Kell] versions of 
molecular codes. These show an improvement of the ionization energies compared 
to Hartree-Fock and Density Functional Theory (DFT). Still, there remain several 
problems with this procedure. 

1. The inclusion of the anti-symmetric Coulomb matrix elements (Generalized 
GW (GGW)) seems to deteroriate the results, while this is essential both for 
extended and finite systems. 

2. The coupling of the single-particle states with two-particle-one-hole and two- 
hole-one-particle states involves only two propagating electrons. This inher- 
ently violates the Pauli exchange with the third electron. 

3. In the GWA, only partial diagonalization in the particle-hole channel is per- 
formed. In doing so, only the particle-hole interaction is screened, while in 
the low density limit this leads to the failure of the GWA. In this limit it 
is important to take into account the particle-hole and particle-particle (pp) 
channels on the same footing. 

4. The best results are obtained by doing calculation on the G W level. This 
means that increasing self-consistency actually makes the results worse, in- 
stead of systematically increasing the accuracy. 

The electron gas in 2D and 3D has also been extensively studied using quantum 
Monte Carlo methods. Provided that the sign problem is under control, these 
accurate methods provide exact results or at least exact constraints on the prop- 
erties of the electron gas that can be used as guidance for the development of 
other approximations or as reference values. More information on the continuous 
formulation of QMC can be found in Ref. [Cep04] and references therein. 




(1.1) 
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Figure 1.1: This figure shows the low momentum behavior of the polarization propagator 
for the uniform electron gas for the Wigner-Seitz radius[DN05] r s = 2. This 
value indicates that the electron gas is in its metallic regime. The discrete 
pole energy versus Q in the TDA is given by the blue diamonds, while the red 
crosses represent values obtained with the RPA. The true plasmon frequency 
should be at -§- = 1 as obtained with the RPA. In the area between the two 
black lines, the polarization propagator has continuous poles. The derivation 
is given in Appendix A. 
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The electronic behavior of solids bears great resemblance to the electron gas prob- 
lem. Consequently, the same methods are applied to the calculation of the selfen- 
ergy and spectral function of solids. The GWA describes the long-range behavior 
of the system in an edequate manner, but lacks the correct short-range interaction. 
A natural extension of the GW formalism therefore includes these interactions by 
using the GW Green's function in T-matrix theory, as such describing multiple 
scattering between two holes or particles[SAK98]. This improves the description 
of the satellite structure of the spectral function, allowing for the occurence of a 
satellite a few electronvolts below the main peak, which is also found in experi- 
ments. This success gives a hint that the simultaneous inclusion of particle-hole 
and particle-particle scattering is mandatory. 

1.1.1b Nuclear Matter 

Nuclear matter (i.e. a homogeneous gas of nucleons) is a system for which the 
particle-particle channel has been studied extensively. The interaction between nuc- 
leons is very repulsive at short distances and the probability of finding two nucleons 
close together is consequently small. As a result, a mean-field description starting 
from the bare nucleon-nucleon potential is not meaningful. However, considering 
the interaction between particles at a more sophisticated level, one creates correl- 
ations that lead to an effective interaction. The dressing of the interaction can be 
attained within the so-called ladder approximation[Gal58]. In this approximation 
only graphs with one hole line are retained and an infinite sum of all ladders is per- 
formed. This approximation corresponds in fact to the low density limit of nuclear 
matter. 

1.1.2 Finite systems 
1.1.2a Molecular physics 

For closed shell atoms, the spectrum obtained in photoionization processes is fairly 
simple. There is a main peak for every orbital that carries the larger part of the spec- 
tral strength, while the rest of the strength is placed in a number of smaller satellite 
lines. For these systems, a single-pole approximation describes the ionization phys- 
ics respectably well. However, in some cases this quasi-particle picture breaks down. 
The breakdown mostly occurs when the occupied levels are energetically close to 
the vacant levels. Whereas for atoms this breakdown of the quasi-particle picture 
is rather exceptional, it is frequently seen in molecules[CDSN07]. This behavior is 
due to the small gap between occupied and unoccupied orbitals in some molecules 
and the large number of occupied orbitals in the valence region. As a result, there 
is mixing of higher excitations with the single-particle states, which implies the 
breakdown of the quasi-particle picture. 
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A successful approach to go beyond the quasi-particle picture is by including two- 
particle-one-hole (2plh) and two-hole-one-particle (2hlp) configurations into the 
eigenstates of the anion or kation[SC78]. This approximation is called the two- 
particle-one-hole Tamm-Dancoff approximation (2plh-TDA) and can be seen as a 
configuration interaction expansion including the 2plh and 2hlp states into excited 
states or as a Green's function method with a selfenergy approximation that ex- 
hibits the right analytical properties of the exact selfenergy. This method leads to 
qualitatively correct ionization properties[CDS80], but takes the ground-state cor- 
relations only partly into account. As a result, outer valence ionization energies are 
systematically too low. Such behavior is cured by refining to an approximation that 
is exact up to third-order of perturbation theory. The extended 2plh-TDA[WS81] 
can be shown to coincide with the third order algebraic diagrammatic construc- 
tion (ADC(3))[SCW83] technique. From the diagrammatical derivation of this 
technique, it is clear that, in the selfenergy of ADC(3), interactions between two 
lines are included up to the Tamm-Dancoff (TDA) level. This method yields very 
satisfactory results for a wide range of small molecular systems[vl\ISC84], but the 
discussion on the plasmon pole shows that the restriction to TDA excitations will 
possibly lead to disastrous results when applied to extended systems. 
To deal with the computational cost of an all-electron calculation, the electron 
ionization and attachment are sometimes treated separately. This leads to the 
so-called non-Dyson ADC (nD-ADC) approximations[STS98]. In this technique 
the ionization and affinity parts are analytically decoupled from the beginning, 
prohibiting the logarithmic divergence of the static part of the selfenergy as shown 
in Ref.[DSC95]. The separation of addition and removal parts of the Green's 
function calculation remedies the violation of particle number that causes this 
problem. To make the method applicable to large bio-molecules, a more drastic 
simplification was introduced[Ort98]. Here, the terms that are unimportant for 
the ionization of closed-shell molecules are discarded. The resulting nondiagonal, 
renormalized second-order (NR2) theory presents a significantly better scaling with 
particle number than ADC(3). 

1.1.2b Nuclei and atoms 

Based on the successes of the techniques using RPA-like interactions and the in- 
dications that coupling to higher excited states have to be included in the theory, 
there have been efforts to unify these qualities in one and the same technique. 
This unification proved to be an extremely challenging problem to solve. Exten- 
sions of the second-order diagram of the selfenergy[RAD92] stay limited to either 
the particle-hole or particle-particle channel. Seeing that both these channels are 
important in some cases[AH71], the inclusion of both channels on the same footing 
was investigated [SVR73]. This resulted in a matrix eigenvalue equation that is very 
similar to the simple RPA equations. These equations were formulated in terms of 
the solutions of the ph and pp RPA. On close inspection, however, there are some 
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severe problems with this approach[RGAD96]. It was shown that the technique 
used to derive the 2plh RPA leads to the inclusion of unphysical solutions. Never- 
theless, the results of the calculations were encouraging and particularly the proton 
spectral strength at low missing energies for 48 Ca nuclei was closer to the experi- 
mental value than the 2plh-TDA result. This was thanks to the extra ground-state 
correlations that induce an extra energy-dependent diagram in the selfenergy that 
is absent in 2plh-TDA. The extra diagram results in a stronger coupling of the 
phonons for states where both the ph and pp channels are important. 

In a recent work[Bar02] a formalism has been introduced that includes the ph and 
pp channel at the level of the RPA. The Faddeev random phase approximation 
(FRPA) starts from an irreducible 2plh/2hlp propagator for the selfenergy, so the 
Pauli principle is fully accounted for. The breakthrough that enables the coupling 
of ph and pp channels on an equal footing is the use of a Faddeev technique[Fad61] 
to resum the diagrams to all orders. The Faddeev route is familiar when solving 
the three-body problem. This results in an improved description of the long range 
correlations that need to be included in the single particle propagator. In this 
approach, the spurious solutions also appear[AG79], but they are easily separated 
from the physical solutions[NBG99] for the three-body problem. In a many-body 
calculation, the separation is less simple but can still be performed [Bar02]. Another 
advantage is that this method is easily extended with different interactions. If one 
wants to go beyond the standard RPA level, this can easily be done by replacing 
the corresponding intermediary excitations[BD03]. The FRPA approach has been 
successfully applied to nuclei[BD02] and atoms[BVND07]. Fragmentation effects 
were also investigated for nuclei. The results obtained with this dressed FRPA are 
in good agreement with the experiments and for some effects it is advisable to use 
the dressed propagators. Overall, the FRPA is a method that is applicable both 
for finite and extended electronic systems because of the RPA screening of the 
interaction. 

In this work we will study the applicability of the FRPA for a set of small molecules 
[DVNB11, DNB11]. In view of the good results obtained in ADC(3) and the 
superiority of the RPA over the TDA, we expect that this method will result in 
a better description of the excitation energies and ground-state energies. We will 
also address the accuracy of the method when the basis set is enlarged[BVND12]. 
To do this we will extrapolate the results for a series of atoms to the basis set limit. 
At every time we will compare our results to the experiment where available, and 
to the current golden standard in molecular and atomic calculations, the CC with 
singles, doubles and perturbative triples (CCSD(T)). 

As an alternative reference method one could also look at the QMC results. It 
is computationally possible to obtain ground and exited state information about 
small nuclei up to 12 C[Wir09] and even in atomic and molecular physics[Cep96] 
QMC methods are becoming more and more established. 



8 



1. Introduction 



1.2 Outline 

In Chapter 2 we will present all the theoretical tools that are needed to construct the 
FRPA. The definition of the single-particle Green's function . The Dyson equation 
will be given in the form that requires the irreducible 2plh/2hlp propagator. After 
this the definition of the polarization propagator and two-particle Green's function 
will be given. These propagators will be used in a suitable approximation called 
the RPA. 

Chapter 3 will deal with the derivation of the FRPA mechanism. For this the six- 
point vertex function has to be reduced from a six-time object to a two-time object. 
The reduction requires that the propagators are two-time quantities and that the 
class of diagrams that is resummed is restricted. The Faddeev procedure increases 
the matrix dimensions by a factor of three. However, by elimination of the spurious 
solutions, the matrix dimension brought back to the original dimension through a 
simple projection. The ADC(3) is encompassed by the FRPA as a limit case which 
can easily be derived. 

The results obtained with the previously derived method will be given in Chapter 4. 
We have applied the FRPA to a series of closed-shell atoms, a set of simple diatomic 
molecules and a schematic model. Through this schematic model we try to give an 
understanding of the problems that occur in the dissociation limit with the FRPA. 
In Chapter 5 we will present the conclusions of this work, together with an outlook 
on possibilities for future research. 



Chapter 2 



Methods 



As the amount of information that needs to be stored for the exact description of 
the wave function of an TV-body system explodes exponentially with the number of 
particles, there has always been a search for other physical quantities that represent 
the physically relevant part of the wave function. Similar to the electronic density 
and the reduced density matrix[Verl2], the single-particle Green's function is one of 
these quantities that describe the physics of an iV-body system in a more compact 
manner. Green's functions form a natural way of describing the evolution of an 
A^-body system. Just like the density matrix, they hold information about the 
ground state and allow for the calculation of the expectation value of any one- 
body operator, as well as the ground-state energy when interactions are limited 
to two-body interactions. At the same time, the poles and amplitudes of the 
Green's function offer access to the N — 1- and A~ + 1-particle systems. For finite 
systems, all data contained by the single-particle Green's functions can be related to 
experiments. The addition process can be probed by elastic scattering experiments 
and the removal process can be witnessed experimentally in (e,2e) experiments. 

2.1 Description of the system 

The time-independent Schrodinger equation can be written as an eigenvalue equa- 
tion 

h\v») = e \*Z), ( 21 ) 

where H is the Hamiltonian of the system and the energy Eq is the eigenvalue of this 
operator for the ground state \^q) of the A^-particle system. The nuclear positions 
change on a much slower timescale than the electrons. One can assume that the 
fields generated by the ions act as an external potential for the electrons. After 
separating the nuclear and electronic wave functions using this Born-Oppenheimer 
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approximation, the electronic Hamiltonian can be written as 

H = ^(a\T\P)aiap + M^a'a* a,*,, ( 22 ) 

where a\ x creates a state with quantum numbers a and a a destroys such a state. 
The quantum numbers for diatomic molecules will typically consist of an index la- 
beling the Hartree-Fock orbitals, a spin projection and the projection of the angular 
momentum on the molecular axis. The matrix elements (a\T\f3) are the matrix 
elements of the kinetic energy operator and the external potential, while the matrix 
elements (a(3\V\"f5) are the expectation values of the anti-symmetrized two-body 
Coulomb interaction. Unlike in nuclear problems, this Hamiltonian fully describes 
the non-relativistic system. There are no three-body forces. Since in electronic 
systems the interaction is completely know, the two-body matrix elements can be 
written in terms of exact 2-electron integrals 

(aftVW) = /dxi y dx 2 4 ( Xl ) ^ (x 2 ) |r ^ r2| 7 ( Xl ) fa (x 2 ) 

- J dxi J dx 2 a (xi) (jyjj (x 2 ) ^ fa (xi) 7 (x 2 ) . 

(2.3) 



The integration over the variables xi and x 2 also implies the summation over the 
discrete degrees of freedom. 

2.2 Single particle Green's function 

To define the single-particle Green's function, we have to switch from the Schrodinger 
picture of quantum mechanics to the Heisenberg picture (Appendix B). In this pic- 
ture the creation and annihilation operators are time dependent 



atjt) = «x,,( ; ' ; ///),/;, «x,,f ; '/// ) (2.4) 



(t) = exp ( l -Ht ) a Q cxp ( -^Ht ) . (2.5) 



The single-particle Green's function is then defined as 



G a At,t') = -^(K\r[a aH (t)4 H (t')}\^), (2.6) 

where l*^) is the exact ground-state wave function from Eq. (2.1). The operator 
T is the time-ordering operator, which introduces a minus sign when t 1 > t 

(*)<(*')] =e{t-t')a aH {t)a) jH {t') - 6{t' - t)a\ H {t')a aH {t). (2.7) 
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The Hamiltonian operators appearing in the exponent of Eq. (2.4) and Eq. (2.5) 
can now be resolved by introducing a complete set of eigenstates for the N + 1 
and N — 1 systems. The Green's function has the form 

x (^\a a \^)(^ +1 \4\<) 

n ^ 

x«| a t|^-i)(^-i| aQ K)} (2.8) 

= G aP {t- If), (2.9) 

which only depends on the time difference t — t'. Note that here we only include 
discrete eigenstates. Taking into account continuum eigenstates would result in the 
replacement of the discrete sum with an integral. The complete set of eigenstates 
for the N + 1 and N — 1 systems have eigenenergies analogous to Eq. (2.1): 

H\^ +1 ) = E^Z +1 ) (2-10) 
H\^-i) = E"- 1 ^). (2.11) 

The first term of Eq. (2.8) is the forward propagating or " particle" part; it describes 
the evolution of an additional particle on top of the TV-particle system with quantum 
numbers (3 created at time t' to a state a which is destroyed at time t > t' . 
The second term is the backward propagating or "hole" part which describes the 
probability of the removal of a particle in state a from the system at time t and 
the return to the ground state at time t' > t by adding a particle in state (3. 
By inserting these complete sets, it becomes clear that the Green's function holds 
information about the N + 1- and N — 1-particle system. This becomes even 
more obvious when we make the Fourier transform to the energy domain where we 
recover the Lehmann representation of the Green's function 

Ga " {E) - \ E-iE^-E^ + i, 



(2.12) 



= G> fj (E) + G<p(E), (2.13) 

in which -q is a small positive number that ensures convergence of the Fourier 
transform. The forward and backward part are labeled in Eq. (2.13). The quantity 
that is probed in (e,2e) experiments is the hole spectral function which is related 
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to the imaginary part of the backward propagating part of the Green's function by 
using the Sokhotski-Plemelj relation 

S<(E) = -^G< a (E) (2.14) 

= £ K*" "Val*^! 2 * - « - ■ (2-15) 

n 

The hole spectral function is a sequence of <5-functions located at the poles of the 
Green's function, which coincide with the transition energies to the N — 1-system. 
At this energy, the spectral function measures the probability of removing a particle 
with quantum numbers a from the ground state while leaving the residual system 
in a state of energy Eq —E. The deviation from the non-interacting single-particle 
picture is best captured by the spectroscopic factor of an excited state n 

z n = EK*"" 1 ! "!*^! 2 - (2.i6) 

a 

This spectroscopic factor expresses the total probability of removing a particle from 
the system and arriving at the excited state n of the iV — 1-body problem. 
The expectation value of any one-body operator can be expressed in function of 
the Green's function. This can be seen from the connection between the density 
matrix and the Green's function. The ground-state expectation value of a single- 
body operator is obtained as 

= ^(aMPWZlciaeW) (2.17) 

a/3 

= ^2(a\0\P)N a p, (2.18) 

afj 

where the density matrix was introduced as 

N afj = (*»\a a a \*»). (2.19) 

By inserting a complete set, the structure of the numerator of the Lehmann rep- 
resentation of the Green's function appears. The single-particle density matrix is 
related to the Green's function by 

n 

- /2d eXpM ^ E — (Eq — En -1 ) — if] (Z21) 
= J ^exp{i V E)Gp a (E), (2.22) 
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where the factor exp (ir/E) is inserted to select the backward going part of the 
propagator when integrating over the complex poles. 

One would expect that the determination of the ground-state energy is beyond 
the scope of the single-particle Green's function, as the expectation value of the 
Hamiltonian involves the evaluation of both one- and two-body operators with re- 
spect to the ground-state wave function. The ground-state energy can be retrieved, 
however, with the Migdal-Galitskii-Koltun sum rule 

< = (9g\H\*») (2.23) 
= V- dE^({a\T\p)+ES^)?sG< a (E), (2.24) 

which greatly extends the applicability of the Green's function theory. Here the 
integral over the energy goes from — oo to the Fermi-level of the system which 
is defined as 

e F = E^-E^- 1 (2.25) 
4 = (2-26) 

The Fermi-level also serves as a separation between particle (p) and hole (h) states. 
The particle states are those with energy larger than or equal to Cp, whereas the 
hole states have energy below or equal to e~j?. 

It is instructive to look at the Green's function for the system without correlations. 
This is the solution of that part of the Hamiltonian H that can be diagonalized 
in a single-particle basis 

H = H + H 1 . (2.27) 

There is a certain degree of freedom for the choice of the operators H and H\. 
One can for instance include the average two-body interaction (mean field) into 
the one-body operator by defining an auxilary potential 

H = T + U (2.28) 

and 

Hi = V — U. (2.29) 

The splitting should be done in a way that the remaining part of the two-body 
interaction Hi is small in some sense, and perturbation theory can be more easily 
applied starting from the eigenstates of H . The operator H Q can be written in its 
eigenbasis as 

H = J2 e »ala a , (2.30) 
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where the e a are the eigenvalues so that 

H \a)=e a \a). (2.31) 

The N-particle ground state of this Hamiltonian is simply the product state of the 
N lowest states en 

N 
i=l 

Inserting this ground-state wave function and Hamiltonian into the definition of the 
Green's function (Eq. (2.12)), one arrives at the non-interacting Green's function 

G (0) 

G%(E) = S a J^^ + ^-^) (2.33) 

p \ h — e a + it] hj — e a — ir\ I 

= G^ f} > (E)+G% < {E). (2.34) 

G^g is diagonal in the eigenbasis of Hq and has a single pole at energy e a with 
transition probability equal to 1. This is also observed in the hole spectral function. 
Large deviations from this kind of behavior in a Green's function indicate a strong 
dissimilarity with a non-interacting system and thus important correlation effects. 



2.3 Polarization propagator and N±2 Green's func- 
tion 

Similar to the definition of the single-particle Green's function in Eq. (2.6), a 
propagator for the evolution of two particles can be defined 



(*ai t(3, t~/i ts) 



(2.35) 



The full G 11 is not needed for the study of excited states. It suffices to take only 
the two-time limit of the two-particle propagator, with the indices corresponding 
to particle-hole (ph) excitations 



lim lim . G^p (t,i? ,tp,t 7 ) 



i 



«|T \aUt)a a Jt)a\ H {t')a- 8H (t')] \^) 



(2.36) 
(2.37) 
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The states a are related to the states a by time-reversal. For an electron in a 
spatial orbital a and with spin projection m s , this means 

|o7m7) = (-l)5+ m =|a,-m s ). (2.38) 

The phase factor in this relation is not unambiguously defined, but it is fixed by this 
definition. Applying the time-reversal relation twice returns the quantum numbers 
to their normal values but with a minus sign 

\a, m s ) = — \a, m s ). (2.39) 

By inserting a complete set of states in analogy to Eq. (2.8), the ph Green's function 
is written as 

G%-i„ S -dt-t) = -^«l4a Q |0«|at a j|0 

-^E^-^cxpQ^ -E»)(t-t') 

x(^|4a a |^)(^|a%-|^) 



X)fl(f-t)expQ«-0(* , -t)) 



h 

x{^\a\a- s \^)(^\ala a \^). (2.40) 

The first term in this propagator still represents a contribution from the ground 
state. It is customary to study the polarization propagator which includes only 
terms that are relevant to the study of excited states and that are not already 
included in the single-particle propagator 

n Q 0-i, 7 i-i(i - *') = G P ap-i tl s-i{t-t') 

+ i«l4a a K)«|ata tf -K). (2.41) 

The non-interacting polarization propagator is again obtained by replacing the full 
Hamiltonian H with Hq and the exact ground state with the non-interacting 
one In the eigenbasis of H this results in 

n S-i, 7 «-i(*-*') = -\s(t-t')e(e a -e + F )e(e- F -ep)5 ai 8ps 



xexp(-^(e /J -e a )(t / -i)) , (2.42) 
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where time-reversed states have been assumed to have the same energy as normal 
states so that the time-reversed notation can be suppressed in the right-hand side. 
The first term describes the simultaneous propagation of a particle a and a hole j3 
at time t' to a particle 7 and a hole S at time t. It is clear that the content of the 
non-interacting polarization propagator can be expressed in terms of single-particle 
propagators of a particle a, making the evolution to a particle 7 and independently 
the propagation of the hole /? to a hole 5. This becomes clear from the formula 

n<g_ 1 . y4 _ 1 (t-t') = -ihG^(t-t')G^(t'-t). (2.43) 

The second term interchanges the time arguments and the roles of the particle 
and hole states. The product of step functions cancels when the time arguments 
have opposite sign. Since the polarization propagator is again only a function of 
a time difference, one can apply a Fourier transformation to get the Lehmann- 
representation as a function of only one energy 



n a/ 3-i i7 5-i(s) = 



(^\ala a \*»)(*»\a\a- 5 \^) 



E 



(*»\a\a s \*Z)(*Z\ata a \*») 



t E+(E^-E^)-i V 



(2.44) 



= K^^-^+K^^iE). (2.45) 

The poles of this function, just like with the single-particle propagator, hold in- 
formation about the excited states of the TV-particle system, while the numerator 
is made up of transition amplitudes. The non-interacting polarization propagator 
in the energy domain is given by 

(0) 1 t?\ XX ( - 4) g «F ~ € fi) 



E-(e a - ep) +ir] 

9{e-F - £ a )0(tp - eg 
E + (ep - e a ) - it] 



(2.46) 



= ^-.^(4 (2-47) 
which consists of a forward and backward propagating term that describe the non- 
interacting ph- and hp-propagation. 

In the same spirit as G ph we can derive the propagator that describes the system 
with a pair of particles removed or added. Again it is not needed to study the full 
two-particle Green's function; the object that holds the interesting information is 
the two-time limit G pp 

G P a P p,-ys(t-t') = lim + l™ Gi^ lS (t,tp,t',ts) (2.48) 

tn— >t+ tx— H' 



1 ■ -,iV 



ap H (t)a aH (t)a\ H (t')al H (t')} K).(2.49) 



2.4. Equation of Motion 
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It becomes clear from the Lehmann representation that we are in fact dealing with 
an object that holds information about the N ± 2-particle system 



G ^ S{E) " \ E-(E»+>-E») + ir, (Z50) 
^ (^\al ( al\^- 2 )(^- 2 \^m 

^ E-{E»-E»-*)-in ■ [t>) 



The non-interacting version of this Green's function is retrieved in the customary 
manner. It can again be expressed as a product of non-interacting single-particle 
Green's functions 



- - (G^(t t')G$(t - - GS(t - f - f )) • (2-52) 



G 



The first term describes the non-interacting propagation of state 7 to a and 8 to 
(3 while the second term describes the exchange equivalent of this, where the two 
initial or final states are exchanged. The propagator in the energy domain shows 
the possible propagation of either two particle states or two hole states 



0(ep - e a )9(ep - epY 



E-(e a + e p ) - ir? 



(2.53) 



2.4 Equation of Motion 



Elegant as they are, exact Green's functions are hard to obtain for realistic systems. 
One has to introduce approximations that describe the correlations in the system 
to keep the procedure computationally tractable. These approximations have to 
be proposed in a way that can easily be expanded to higher orders. One such 
approach is the equation of motion approach for the Green's function. It starts 
from a Hamiltonian written in the form of Eq. (2.27). It is possible to write the 
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time derivative of the single-particle propagator as 



(2.54) 



<dt 



6{t-t')a aH {t)a\ H {t') 
-e{t'-t)a\ H {t')a aH {t)] K) (2.55) 
= S(t tf)(^\a aH (t)4jt') + a\ H {t')a aH {t)\^) 



Of, 



\K) 



(2.56) 



5(t-t)6 af) + {*Z\T 



\K) 



(2.57) 



where the fundamental anti-commutation relation for the fermionic operators was 
used. It suffices to know the time evolution of the annihilation operator a aH (t) to 
determine the equation of motion for the single-particle Green's function. In the 
Heisenberg picture we find 



ih—a aH (t) 



= K H (t),H] (2.58) 

= e a a aH (t) - ^2(a\U\S)as H (t) 

s 

+ lY,( a6 \ V ^04 H (tMtK H (t)- (2-59) 

<5,7,C 

Substituting this in Eq. (2.57) we get 

G a f}{t — t') = 5{t-t')5 a p-Y,{(AU\S)Gs^t-t l ) 

s 

+ lY,( aS \ V MKi S e(t,t,t + ,t'). (2.60) 

<5,7,C 

This equation relates the single-particle propagator to the 2-particle propagator. 
Following this path, it is possible to express every iV-particle propagator in function 
of the N+ 1-particle propagator, building a hierarchy of more complicated objects. 
Writing the equation of motion for the non-interacting single-particle propagator 
yields an equation similar to Eq. (2.60), but without the interaction terms 



lh dt £c 



lh dt tc 



G<$(t-t') = 5(t-t')S 



(2.61) 



2.4. Equation of Motion 
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This defines the inverse of the non-interacting propagator as 



G«8~\t-t') = 5(t-t')8 



■n 9 

in— - e c 
ot 



(2.62) 



which is exactly the first factor in Eq. (2.60). Bringing this factor to the right-hand 
side of the equation leads to 

G aP {t-H) = G^(t-t')- JdtiY.GMit-hMUlQGspih-t) 

7,<5 

+\ I dti J2 G^(t-h)(iS\V\ 7 

J <5,7,C,<. 

xG^fi^.t+t'). (2.63) 

One can now split off the reducible parts from the two-particle Green's function as 
a product of exact Green's functions and replace the interaction with the so-called 
four-point vertex function r 4 ~ p *, which serves as an effective interaction between 
dressed particles in the medium. It has four time arguments, even though the last 
term in Eq. (2.63) only has two time arguments. The connection to the external 
times is also handled by interacting single-particle Green's functions, which results 
in the following expression for G 11 

G^ i74 (t,t,t+ ,t') = iHG a7 {t-t+)G l3 6(t-t , )-ihG a s(t-t')G 01 (t-t + ) 

+{ihf ( dt e f dt c ( dt, f dt G <*e(t-t e ) 

x G K (t - t c ) «|r 4 - pt (t e , t c , t„ t e )\rfi) 
xG^-^Gesite-t'). (2.64) 

Inserting this equation in Eq. (2.63), one can define the irreducible selfenergy E* 

£*/,(*, i 7 ) = -5(t-t'){a\U\p) -ih6(t-t')J2(ai\V\pS)G S7 (t-t+) 

7,(5 

+ (ri<f\ j dt ( / dt, J dt e Y, (uj\V\5e)Gs C (t - t ( ) 

7,<5,e,£,r/,f9 

xG en {t-t 11 ){(:i 1 \Y^ t {t t 11 ,te,t')\ei3)Ge 1 {te-t), (2.65) 
so that the full Green's function is the solution of the Dyson equation 

G afj (t-t>) = G^(t-t') + j dt, j dt 2 Y G ^(t-ti)^ 5 (h,t 2 ) 

7,5 

xG 5 p(t 2 -t'). (2.66) 
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This is a self-consistency problem. The full Green's function appears in the left- 
and right-hand side of the equation, as well as in the irreducible selfenergy. The 
correlated Green's function is equal to the non-interacting Green's function and 
a part that includes all interactions. Making an approximation to the propagator 
implies finding an approximate selfenergy. 

The irreducible selfenergy S* acts as a potential for the Green's function. This can 
be made clear by going to the energy domain, where the convolution time integrals 
of Eq. (2.66) turn into simple products in the energy representation: 

G aP (E) = G^{E) + -£G^(E)Ti* s {E)G S p(E). (2.67) 

7,(5 

Multiplying by a factor E — e n , with e„ a pole of G in the removal part of the 
Lehmann representation, which naturally does not coincide with a pole of G^ ' or 
S* and then taking the limit 



lim (E — e„ 



G a p(E) = G { a °>(E) + ]T G^(E)T^ S (E)G SP {E) 



7,6 



(2.68) 



and dividing by a factor |al|\E'^ r 1 ), one arrives at an eigenvalue equation for 
the amplitudes of the Green's function 

(eVolO = ^GW(e„)S; 5 (e„)(^- 1 |a 5 |^). (2.69) 

7(5 

This was derived working in the basis that diagonalizes H , but Eq. (2.69) is valid 
in any basis. 

Special attention should go to the normalization of the left-hand side of Eq. (2.69). 
One can prove[DN05] that the states are normalized like 

<W-^f(e„) W|at|^-i) = L (2?0) 

The square of the norm of these amplitudes constitutes the normalized spectro- 
scopic factors that are commonly determined in experiments employing (e,2e) or 
(e.e'p) reactions. 

Now one can make the transition to the coordinate representation a = rm s 
<*rW s lO = E / dr i/ dr 2 GW sri (e n ) 

mi, 7712 

x^ imir2m2 (e„)(^-Vr 2TO2 |^). (2.71) 



2.5. Diagrammatic language 
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Using the coordinate representation of Eq. (2.61) 

^2 J ^{rm s \e n -H Q \r x m 1 )G { ^ mir , m , s {e n ) = <5 msI< 5(r - r'), 

mi 

(2.72) 

and writing the same equation for the amplitudes ( v f'^ r_1 \a rmg I^q) 
/ dr 1 (rm ;s |e„ - i? |r 1 TO 1 )(*^- 1 |a rimi |*^) 

mi 

= ^ + ~ U(r) ] ^n-^rmAK) (2-73) 

generates a Schrodinger-like eigenvalue equation 

-^(^KmJ^+E / dr 1 S^ srimi (e„)(^- 1 |a rimi |^) 

mi 

= £n (*r i |«rmjO(2.74) 

where S* represents the irreducible selfenergy without the external potential. From 
this equation it is clear that the irreducible selfenergy acts as an effective poten- 
tial for a single particle and takes into account the interactions of the particle 
with the medium. The same equation can be derived for the forward propagating 
amplitudes, with the understanding that these states with e„ > ep~ correspond to 
scattering states. 

2.5 Diagrammatic language 

One of the great advantages of Green's function techniques is the existence of a 
graphical language to express complex equations. Every element of the Green's 
function theory can be represented by a Feynman graph[MP92, DN05]. Table 2.1 
presents the propagator and interactions in the time domain that are used in Sec- 
tion 2.2. In each of these diagrams, both the forward and backward term are 
represented in one diagram, i.e. there is no time ordering. The full propagator is 
presented by two lines, while the non-interacting propagator is represented by a 
single line. 

The same procedure can be repeated in the energy domain. Every propagator 
describes the evolution of a state with energy E, both forward and backward 
propagation are captured by the same diagram. 

The Dyson equation Eq. (2.67) is depicted in Figure 2.1. From this diagrammatical 
representation it can easily be seen that the Dyson equation represents an infinite 
resummation of diagrams: this infinite repetition acts as a renormalization for the 
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Diagram Time domain Energy domain 



G ( >-t') G^JE) 



G a p{t — t') G a p{E) 
(a\U\p) (a\U\p) 



(ap\Vh8) (af3\V\ 7 S) 



Table 2.1: Diagrammatic translation of the propagators in the time domain and energy 
domain 





Figure 2.1: Diagrammatical representation of the Dyson equation in the energy domain 



2.6. Conserving approximations for the selfenergy 
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(b) 

Figure 2.2: Two different ways to write the selfenergy E* as a functional of the four-point 
vertex function r 4_pt 

interaction. While the Dyson equation is in principle exact, the determination 
of the complete selfenergy is just as hard as that of the full propagator. In most 
approximation schemes for the selfenergy, the number of diagrams that is resummed 
is restricted to a certain class of physically relevant diagrams. 

2.6 Conserving approximations for the selfenergy 

Eq. (2.65) has a diagrammatical form as displayed in Figure 2.2a. There is, how- 
ever, another choice possible in the definition of the four-point vertex, displayed in 
Figure 2.2b, or as an equation 

53^(4,0 = -5{t-t'){a\U\P)-ih5(t-t')Y J {ui\V\P8)Gs 1 (t-t+) 

7<5 




xG v5 (t v - t')G Se (t e - t')(6e\V\p 7 )G Cy (t - t c ). (2.75) 

There is no reason why the two different definitions should lead to a different 
outcome. In fact, it has been shown[BK61, Bay62] that the equivalence of these 
two diagrams has important physical implications. If these are equivalent and the 
four-point vertex is symmetric 

^a/3 P -yS ~ ^p a P 5f (2-76) 

then the particle number, total momentum, total angular momentum and energy 
are all conserved quantities in this approximaton. This is a very strong statement 
and approximations that have this property are to be preferred above others. 



24 



2. Methods 



To enforce this symmetry between the two ways to write the selfenergy, it is best 
to study the equation of motion for G 1 ^ 7<5 (i, i, t + , t'), but now with the time 
derivative with respect to t' 

d 

- ih ^rr,Gai3 ijS (t,t,t + ,t') = (S a ^8 fjS - S aS 5 fj7 ) 8(t - t') 



dt 

d 
dt 



+(*o\T 



(2.77) 

using the equation of motion for the creation operator as H (t'), we arrive at 
(-ih-^ - e^j G I J /3 lS (t,t,t + ,t') = (S a7 5ps - S a sSpj)6(t-t') 

-£G^ i7e (M,i + ,i'KW> 

e 

-^E g SK«c,(*-«')«m*»?> 

eC>7 

(2.78) 

which is the definition of a next Green's function in the hierarchy G 2plh 

(2.79) 

This Green's function describes the simultaneous motion of two particles and one 
hole in the medium. Just like the single-particle Green's function, it depends on the 
difference of two times only. The first factor in the right-hand side of Eq. (2.78) 
represents a non-interacting propagator G^°\ The right-hand side can be rewritten 
in terms of the second definition of the selfenergy Eq. (2.75). Using the Dyson 
equation Eq. (2.66) this results in the following expression for G*Jp lS (t 1 t,t + ,t') 



G I J 0nS (t,t,t + ,t') = \ ]T fd^Gl^it-t'XeClV^Gesih-t') 
~\ S / dtl J ^dhG^^t^t+M) 



e,C,7J,$,t,K ' 

x G £ - 1 (t ! - t 2 ) Gg iflt (t 2 , *3 , t 3 , t 3 ) (7,0 \V | Ki) 
xG Kg {t 3 -t') 

+ J2 VtrrSfc - SaeSfr) G eS {t - t'). (2.80) 



2. 7. Hartree-Fock 
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Figure 2.3: Diagrammatical representation of the selfenergy from Eq. (2.81) that fulfills 
the Baym-Kadanoff requirements 



Substituting this in Eq. (2.63), one gets a new definition of the selfenergy in terms 
of 2plh excitations[ES69a, ES69b], which automatically fulfills the Baym-Kadanoff 
requirements 

Kp(t-t') = -5(t-t'){a\U\p)+5(t-t')Y d {^\V\p8)G l5 {t-t + ) 

+ Y, (ailVlSejRySeXrteit-t'mvlVWe). (2.81) 

This is graphically represented in Figure 2.3. The irreducible 2plh/2hlp propagator 
R represents the propagation of two-particle-one-hole states and two-hole-one- 
particle states without intermediary annihilation of a particle or hole state 

iWcc(*-0 - Gl% h S< (t-t')-J2 J dti J dt a G I J PntV {t,t,t+,t 1 ) 

7],6 

xG^ii- *2)G^ eC (t 2 ,t'-, t',t') (2.82) 



2.7 Hartree-Fock 

One can in principle choose any approximation for the irreducible 2plh/2hlp 
propagator in Eq. (2.81), but solving for G becomes too difficult very quickly. 
The easiest choice is just to set R equal to zero. In that case, only the first two 
diagrams in Figure 2.3 are included in the selfenergy. This approximation is called 
the Hartree-Fock (HF) approximation. By construction, this method satisfies the 
Baym-Kadanoff requirements. 

E*f(t-f) = -S(t-t')(a\U\/3)-ih6(t-t')Y(^\V\l3S)Gf 7 F (t-t+) 

(2.83) 

This equations explain the name self-consistent field equation, the Hartree-Fock 
quasiparticles move in the field generated by all the other quasiparticles. The 
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Figure 2.4: Diagrammatical representation of the Bethe-Salpeter equation for the po- 
larization propagator 



selfenergy that is generated in this approximation is an energy-independent one- 
body operator. The energy independence of the selfenergy implies that this Hartree- 
Fock prescription is closely related to an independent particle picture. One can 
construct a basis in which the propagator is diagonal. 

The choice of Hq in Eq. (2.27) is fixed only up to an auxilary one-body operator 
U. The choice of U can in principle be completely arbitrary, but it is best to 
make a choice that already reflects the physics of the many-body system. Since 
the HF Hamiltonian is a one-body operator only, it can be chosen as Hq. The HF 
propagator can be seen as a suitable candidate for G^ - 1 . In the rest of this work, 
whenever is used, it is assumed to be the HF Green's function. 



2.8 Random Phase Approximation 

Finding a suitable approximation for the polarization propagator II starts with the 
Bethe-Salpeter equation[SB51, Nak69], which is like a Dyson equation for the 
polarization propagator. This diagram is shown in Figure 2.4. In principle the 
interaction K^' 1 ^ can be energy-dependent, but in what follows we use a simpler 
approximation. The logical first step in this approach is to take K^ ph ^ equal 
to the anti-symmetric two-body interaction V . The polarization propagator that 
is reached in this way is called the Random Phase Approximation (RPA)[BP53] 
polarization propagator. It represents an infinite resummation of all ring or bubble 
diagrams. The argument that the diagrams that are not present in this sum have 
a small contribution due to a randomly changing phase factor and the forthcoming 
cancellation of terms is the historical reason for the name. 



2.8. Random Phase Approximation 
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Figure 2.5: Diagrammatical representation of the Bethe-Salpeter equation for the RPA 
polarization propagator 



The self-consistent equation in the energy domain is given by (see Figure 2.5) 

+ E n 2-,a^ (£) K>l/3e>n£ P V 1 ( E ) ( Z84 ) 

and represents both forward and backward propagating rings. In this approximation, 
it is possible that multiple particle-hole excitations are present at the same time. 
It is convenient to write the amplitudes for the propagator as 

4f )n - <^kM<r (2.85) 

y% k)n = Mato^)*, (2.86) 
so that the RPA propagator becomes 

U RPA rP N \ - ^afi V^y5 > \ - U>og J ^rj , 9 R7 x 

n^-^OE) - >") +Mi r + >") in ' (Z87) 

where the pole energy is written as 

eg*> - £^ - <. (2.88) 

To calculate the poles and amplitudes one has to resolve to a technique introduced 
in Eq. (2.68). This results in a non-hermitian eigenvalue equation 

A B\( X^ \ {ph) (I \ / XW* \ , . 

v b* A* J \ y(p h > ) ~ » ^ o -l ) \ y(p h > J ' [ ' 

with 



Aph,p'h' — fipp'fihh' (£p e h ) + (ph'\V\hp') (2.90) 
Bphjh' = (pp^M'), (2.91) 
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where the indices p and h obviously indicate a particle and a hole state. The 
solutions are normalized according to 

)n r^t )n ' - Y, ( y% h)n )*y P t )n ' = <w. (2.92) 

p,h p,h 

In the same way there is also the closure relation 

E 4 p h h)n ( x X' n y E y% h)n (y { fv n r - v<w. (2.93) 

n n 

This approach differs from the Tamm Dancoff approximation (TDA) in incorpor- 
ating correlation also into the ground state, whereas the TDA takes into account 
correlations only for the excited states. The difference in terms of diagrams is 
represented by Figure 2.6. The TDA propagator does not include backward going 
terms, so there is only one excitation present at any time. In the TDA the matrix 
elements of B are put equal to zero so that the eigenvalue equation reduces to 
a Hermitian problem. This has the advantage that it eliminates the possibility 
of complex eigenvalues when the Hartree-Fock approach reaches instability. The 
result of putting B equal to zero is that there are no backward going amplitudes 
y(ph) i w hich means that these amplitudes can be seen as a measure for the ground- 
state correlations. Due to the lack of backward going terms, the normalization and 
closure in the TDA simplify to 

Y^X^yxW"' = 5 nn , (2.94) 

p,h 

and 

J2 x X h)n ( x X' n r - <w<w- (2.95) 

n 

There is another interpretation of the RPA formalism. These equations can be 
derived by using the equation of motion (EOM) method of Rowe[ROW68]. An 
excitation operator Q( ph ^ has to be defined 

Q 0*)»t = E^^-E^Xap, (2-96) 

p,h p.h 

so that it both creates and destroys a ph pair. Assuming that the hermitian trans- 
posed operator destroys the ground state |5'( p ' 1 '), one can carry out a minimization 
of the energy with respect to the two amplitudes X^ ph > n and y(P h > n . This results 
in the two equations 

(¥ ph )\[ala h , [fr,Q(P h > n t]]|tf(P fc )) = e 0>ft)<*Cpft)|[ a t ahj g(ph)»t]]|^(pfc)) 

(2-97) 

(* (ph) |[4a p , [fr,Q(" h > n t]]|tf(P fc )) = e^tf^l^ap.Q^"*]]!*^). 

(2.98) 



2.8. Random Phase Approximation 
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(b) 

Figure 2.6: The difference in polarization propagators in the RPA and TDA 
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Figure 2.7: The ladder diagrams included in pp RPA. 



When one then again assumes that the ground state is not very different from 
the HF ground state, the double commutators in the previous equations are much 
easier to evaluate. The equations thus derived are exactly the same as Eq. (2.89). 
Replacing the RPA ground state with the non-interacting ground state is called the 
quasi-boson approximation[BET61]. The particle hole operators are considered to 
be perfect boson operators, which is a violation of the Pauli principle, as in fact 
they are built up out of fermion creation and annihilation operators. 
The transition to TDA can also be easily understood from this approach. The 
same procedure can be applied but with the operator 

p , h 

which involves only the creation of a ph pair, disconnecting the ph and hp excita- 
tions. 

The prescription in Eq. (2.96) allows for a straightforward extension of this RPA 
to pp and hh excitations (see Figure 2.7). The excitation operator becomes 

Pl,P2 h 1 ,h 2 

Following the same procedure as for the ph excitations, the result is a non-hermitian 
eigenvalue equation 

/ A B \ ( X^)m \ _ / ! W X ( PP )rn X 

I ijt c J \ y&p)" 1 J m \ o -l J \ y(pp) m J ( ' 

for the eigenvalues associated with the N + 2-particle system, labeled by a plus 
sign and 

( A b \ ( yd*)* \ _ p) _ f i o \( \ 



2.8. Random Phase Approximation 
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for the eigenvalues associated with the N — 2-particle system, labeled by a minus 
sign. The matrices involved are the same in both equations 

) + (piP2\V\ P ' lP ' 2 ) (2.103) 
B PlP2M h 2 = -<PiP2|V|/»i/i2) (2.104) 
Ch 1 h 2 ,h' 1 h' 2 = -8h 1 h' 1 Sh 2 h 2 (eh 1 + eh 2 ) + (/iiM^K^)- (2.105) 

The normalization and closure are very similar to the ph RPA: 

E i-K%l m r x ^l m ' E iy& m )*y& m ' = ^ (2.106) 

E («")*«"'- E (^£ B )^ffi B ' = ~<W (2-107) 

Pi>P2 hi,h 2 

and 

E4? ro (^f )ro )*-E(#")^S p) " = ^(<w^-<ww> 

n m 

(2.108) 

where a, /?, 7 and <5 are generic indices that label two particles or two holes. The 
transition to TDA is also made by setting the backward going amplitudes to zero. 



Chapter 3 



Faddeev Random Phase 
Approximation 



In this chapter we will demonstrate how a selfenergy of the type Eq. (2.81) can be 
built, so that it includes interactions between each pair of states that are at the 
RPA level of theory. This approach will have to go beyond the naive approach that 
sums the diagrams where the intermediary propagator has been replaced by an RPA 
propagator. Such a procedure would introduce unphysical solutions of the Dyson 
equation due to the third diagram in Figure 3.1. This diagram is included in both 
the first and the second diagram and must be subtracted to avoid double counting 
of diagrams. As a consequence of the minus sign in front of the third diagram, there 
will be meaningless solutions to the Dyson equation and the normalization of the 
amplitudes will become problematic due to the erroneous behavior of the derivative 
of the selfenergy. This naive resummation also lacks the exchange diagram that 
lets the first particle line interact with the hole line, which is a violation of the 
Pauli principle. There is no reason why one line should be privileged. 

To cure this behavior, it is necessary to perform a partitioning of the interaction 
into three channels using the Faddeev technique[Fad61]. This allows the pp and 
ph channels to be mixed, while not making a distinction between any of the three 
lines. As a result, a better treatment of Pauli correlations is obtained. In contrast 
to the naive resummation, there is no need to subtract double counting diagrams. 
The Faddeev procedure takes care of the interaction between two lines separately, 
which results in a more elegant simultaneous inclusion of pp and ph channels. 
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Figure 3.1: Naive resummation of RPA interactions between pairs of lines. The third 
diagram introduces spurious solutions in the Lehmann representation of the 
propagator. 



3.1 Bethe-Salpeter equation for the irreducible 2plh/2hlp 
propagator 

It is possible to write a Bethe-Salpeter equation for the irreducible 2plh/2hlp 
propagator[Win72, ES69a, ES69b] very similar to the one for the polarization and 
two-particle propagator. The equation in the time domain 

Ra/3j,Se(;(tl,t2, £4, t 5 , < 6 ) = G a S {t 1 — <4)G^ e (t 2 — t 5 )G(y (t 6 — t 3 ) 



describes the 2plh propagation as the dressed propagation of the 3 states sep- 
arately, including exchange and capturing all interactions between the lines in an 
interaction vertex K . This definition of the interaction vertex can be treated us- 
ing a Faddeev procedure[Fad61], including a separate interaction vertex for every 
channel, as is illustrated in Figure 3.2 

KaPi,5eZ (il)*2)*3)*4>*5)t6) = K^ l ^(t2,t 3 ,t 5 ,t e )G a g(ti—t4) 



As can be seen in Eq. (3.1), the full six-time-dependence is needed to write down 
this Bethe-Salpeter equation. The self-consistency equation can be closed with 
a four-time reduction of the vertex function, but fewer time arguments are not 



—G ae (ti — ti)Gj3s{t2 — t 5 )G^(t 6 — t 3 ) 





(3.2) 



3.1. Bethe-Salpeter equation for the irreducible 2plh/2hlp propagator 
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Figure 3.2: Bethe-Salpeter equation for the six-time irreducible 2plh/2hlp propagator. 
The interaction is split up in separate channels. 
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possible due to the full time dependence of the interaction vertex K . There has 
been no approximation in writing down Eq. (3.1), so if all objects retain their full 
time-dependence, this equation solves the 2plh motion exactly. The K^ pp ^ and 
j((ph) k erne | s describe the irreducible interaction of only two lines in the diagram, 
while the third line propagates freely. It should be noted that K^ pp ^ and K^ ph ^ 
depend on four times. The K^ 2plh ^ kernel describes the interactions between the 
three lines that cannot be separated as a succession of pp and ph interactions. 
In what follows, we will neglect these effective three-body interactions since they 
do not appear at zeroth order in the problems we are addressing, and will be less 
important at higher orders. There are no double counting diagrams so it is not 
necessary to subtract any terms in the equation. 



3.2 Faddeev procedure in the time domain 

It is customary[GI83] to define three Faddeev components i?W that will be used 
to solve Eq. (3.1). The superscript (i) stands for the number of the line which 
propagates without interacting with the other lines, while the superscripts (j) 
and (k) are the lines that are involved in forming the intermediary interaction 
vertex. The triplet (i,j,k) is always a cyclic permutation of (1,2,3). The Faddeev 
components are defined as 

Raly,6etfa>-h,ta,U,tB,t6) = ^ J dt' 2 J dt' 3 J dt'^G f^{t 2 - t' 2 ) 

X Ram,5eC, (^1 > ^2 > *3 ' ^5> *e) (3-3) 

ri,8,L,K 

xRt,/3K,8e({ti, tz,t 3 , ti, t 5j te) (3.4) 

X^K7,5eC(*l)*2! *3, *4 5 h, t&)- (3-5) 



3.2. Faddeev procedure in the time domain 
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Figure 3.3: Every Faddeev component ends in an interaction which lets line (i) go 
freely. Everything that happens before this interaction is captured in the free 
propagation and the other two Faddeev components. 



In the current case where the 2plh kernel K^ 2plh > is not taken into account, the 
full irreducible 2plh/2hlp propagator is given by 

RaPy,5e(;(tli ^2) *3> ^4) *6j — G a $ (ti — to)G p £ (t 2 — tfrjG^ (tg — t 3 ) 

—G ae (ti — te)Gpg(t2 — t5)Gjc{t6 — *3j 

+ X] ^a^7,iec(*l>*2>*3)*4)*5>*6)- (3-6) 
i=l,2,3 

The meaning of the upper index (i) is clarified by the equations 

fl$ 7 ,j e c(*l.*2>*3>*4,*5,*8) = j dtj j dt& j dt9 j dh ° j dhl j ^ 

G arl {ti — tj)Gpg{t2 — t$)Gs L (tg — t 3 ) 

Xr r ) e,,,KA At (*7, ts, td, ho, tiutu) 

X (G K s(tio — t4,)G\ e (tn — t 5 )G^(t 6 — t 12 ) 

— G\s(tio — t4)G Ke (tu — tc,)G^(tQ — tyi) 

+^i*A M ,<5eC^ 10 ' *H>*12i *4j h, te) 



12 



Every Faddeev component i?W ends with an interaction block while everything 
that happens before is captured in the free propagators and the other two compon- 
ents RP' and R^ k \ This is illustrated in Figure 3.3. The interaction blocks can 
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be derived from the exact Bethe-Salpeter equations for the pp and ph propagator 
(Figure 2.4). This Bethe-Salpeter equation can also be written down for the pp 
and ph four-index interaction blocks. The screened interaction will be equal to the 
sum of the lowest-order interaction and higher-order diagrams that are generated 
from a self-consistency procedure. In this way the interaction blocks are derived 
for both channels 



J dt 5 J dt 6 J dt 7 J dt 8 r2J ) eC (ti,t2,t5,*6) 



r{ P h) 



xG er ,{t 5 -t 7 )G ( e(ts-t 6 )K^> s (t 7 ,t 8 ,t 3 ,U) (3.8) 

r^ s (ti,t 2 ,t 3 ,u) = K%%(t u t 2 ,t 3 ,t 4 ) 

+ J2 / d *s / dt 6 f dt 7 f dt 8 r^ ) eC (ti,t 2 ,t5,t6) 

e,C,ri,0 J J J 

y-G cv {t 5 -t 7 )G C e{t & -t & )K^ e %{t 7 MMM)- (3.9) 

The interaction blocks rW are six-index interactions composed of the four-index 
interaction blocks together with the inverse of a single-particle propagator. Two 
of the rW involve the ph interaction, while one is composed of a pp interaction 

= T^ihMMMG-^ih-U) (3.11) 

ri 3 j 7i ^(ti,*2,t3,i4,i5,i6) = r^ ) (Se (ti,t 2 ,t4,i5)G 7C 1 (t 6 - 1 3 ). (3.12) 

Eq. (3.7) still depends on six times and involves integrals over six more times. 
Solving this kind of equations is beyond the reach of present-day computers and 
needs to be avoided at all cost. It should be feasible to circumvent the time- 
dependence of the Bethe-Salpeter equation and to reduce the time-dependence of 
R to only two times. The standard method to reduce the time-dependence of this 
vertex function is by integrating out the excess times. Obviously, this does not 
offer a proper solution. One part of the solution will be to include a polarization 
propagator and two-particle propagator that only depend on two times. The other 
part of the solution is to restrict the class of diagrams that are included to a smaller 
set. 



3.3 Faddeev procedure in the energy domain 

The first step is the reduction of the polarization propagator and two-particle 
propagator to their RPA form, by replacing the kernels K^ ph ^ and K^ pp ^ by their 



3.3. Faddeev procedure in the energy domain 
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lowest-order approximation, i.e. the bare interaction. In the energy domain, the 
ph and pp vertices, like the propagators, will depend on only one energy. They are 
the solution of 



-,{ph)RPA 
a/3, 7(5 



(E) 



VW 7 + 2^ 1 a/3 



(ph)RPA 



J ^ G "' {E + E ^ G ce( E ^e, (3.13) 



-,(pp)RPA 
a/3,-fS 



(E) 



e£,ri,6 



x / ^G^E + E^i-E^V^s. (3.14) 

The vertices can be derived from the polarization propagator and two-particle 
propagator. As a result, their Lehmann representation has the same poles as the 
RPA propagators, but with different amplitudes 



-,{ph)RPA 
a/3, 7<5 
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j,(ph)n _ V^ T/ y{ph)n 
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and 



a/3,7<5 



(3.19) 
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The intermediary single-particle propagators are non-interacting. This is a reflection 
of the quasi-boson approximation in the RPA. In principle dressed propagators[BD02, 
Bar02] should be used in the Bethe-Salpeter equations. This would result in differ- 
ent equations for the polarization propagator and two-particle propagator, which 
are called the dressed RPA (DRPA) equations[RGB+93, Bar02]. The inclusion of 
these effects involves huge matrix sizes as the size of the matrices cubes at every 
iteration. Therefore, in this work there will be no dressed propagators in the Bethe- 
Salpeter equations. Since the polarization propagator and two-particle propagator 
are only calculated with Hartree-Fock propagators, the quality of the vertex func- 
tions would be such that it has no advantage to use dressed propagators elsewhere. 
From now on the non-interacting propagators will be used for the propagation of 
particles and holes between interaction vertices and vertex functions. This guaran- 
tees that at every moment the same level of theory is maintained. In what follows 
it will also be of some technical importance (see Eq. (3.36)) that the poles of the 
intermediary propagation are the same energies as those used in the RPA matrices. 

These spectral representations of the vertex functions will be resummed in the 
irreducible 2plh/2hlp propagator. As a consequence of the reduction to one energy 
in the vertex, the irreducible 2plh/2hlp propagator only depends on three energies 



R^SedE^E^Es) = G^iE^iE^i-Es) G^(Ei)G^(E 2 )G^(-E 3 ) 

+ £ G^(E 2 )G^(-E 3 )V VK ^ 





x / — . 4 ^7,^(^4, Ei+ E 2 - E 4 , E 3 ). 




(3.23) 



For the Faddeev components this implies the same decrease in complexity, to only 
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R^^EuEs) (3.24) 

x j^ i {G^{E 1 )G^XE i )G^{E 2 + E 3 -E i ) 
-G^(E 1 )G { °\E i )G^(E 2 +E 3 - E A ) 

+RaiK,6ec(El, E4,, E 2 + E 3 - E 4 ) 

+Ci, & ((£i^4, E 2 +E 3 - E 4 j) (3.25) 

x f^(G^(E 4 )G^+E 2 -E 4 )G^(-E 3 ) 

-G^>(E i )G^(E 1 + E 2 - E 4 )G^(-E 3 ) 
+- R [k 7 ,5£c( £; 4, Ex + E 2 - E 4 , E 3 ) 
+R[% MC (E^E 1+ E 2 -E 4 ,E 3 ). (3.26) 

While representing a serious decrease in the complexity of the equations, the in- 
tegral over the single energy still forms a mathematical problem that needs to be 
avoided. This means a reduction to a single-energy object is mandatory. 
The reduction to a single-energy irreducible 2plh/2hlp propagator is accomplished 
by restricting the class of diagrams that is resummed in the Faddeev procedure. 
In the discussion so far, the 2plh and 2hlp parts of the irreducible 2plh/2hlp 
propagator R have still been linked. The propagators in between phonons have 
both forward and backward propagating parts, potentially connecting a 2plh state 
with a 2hlp state. This possibility will now be excluded by restricting the number 
of phonons that can propagate to one. There will still be an infinite resummation of 
all possible combinations of ph and pp RPA phonons (including backward propaga- 
tion), but there is only one phonon propagating at every step in time. In between 
the phonon propagation, the three states will propagate forward in time. The omit- 
ted diagrams will only contribute at higher orders and can safely be discarded. This 
procedure unlinks the 2plh and 2hlp spaces and also allows to calculate the irre- 
ducible propagator for the two spaces separately. The spaces will still be connected 
indirectly through the mixing with single-particle space. Even though the 2hlp is 
more directly linked to experimental results, it is important that R is calculated for 
both spaces to include the indirect infuence that both spaces have on each other 
through the coupling with the single-particle states. The derivation of R in the 



one energy integral 

<U C (£l,£2,£ 3 ) = 
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3. Faddeev Random Phase Approximation 




Figure 3.4: An example of a diagram that is not included in the summation of the FRPA. 

This diagram links the 2plh space with the 2hlp space, which is prohibited 
in the present approach. 



2plh and 2hlp space will be very similar, so in what follows the discussion will be 
restricted to the 2plh space. 

Limiting the number of phonons to one at every time step is analogous to replacing 
the factors of three propagators involving forward and backward parts with one 
propagator depending only on one energy. Since from now on the direction of 
propagation is fixed, the indices will indicate whether a state is a particle or a hole 
state. The replacement of the propagators is performed by the substitution 

C;^)^^)^^^) G*> hMh ,(E), (3.27) 

where the > indicates that only the forward propagating part is retained. The 
propagator that describes the evolution of the 2plh state is given by 



r< (0)> frp\ _ fipiP' 1 dp 2 p! 2 Sh 1 h' 1 . . 



Using this propagator for the intermediary evolution of the states results in a new 
equation for the Faddeev components, depending only on one energy argument: 

pW it?\ _ (~i(o)> (E)T^ (E) 

p 1 p 2 h 1 ,p' 1 p'r ! h' 1 \ > / j PlP2h-L,P3P&fl2 \ i P3Pih 2 ,p 5 p 6 h 3 \ I 

P3,P4,P5,P6,h 2 ,h 3 

x (G (0)> h , , hl (E)~G {0)> , , , h AE) 

\ P5P6'l3,PiP2"'l V ; P5P6'l3,P2Pl' l l V ' 

+R pLh 3 , P ' lP ' 2 h' 1 ( E ) + R pLh 3 ,p' lP ' 2 h[ (£)) ■ ( 3 - 29 ) 



3.4 Derivation of the eigenvalue equation 

Fixing the propagation direction of the motion between the phonons has repercus- 
sions on the interaction boxes T^ l \ There are terms that make the different 
from the norma I r^P) and vertices. The real structure of these is needed 
to obtain an eigenvalue equation that is formulated in function of the RPA quant- 
ities and two-particle interaction matrix elements. As an example the diagram 
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Figure 3.5: Diagrams that illustrate the fixing of the time direction. The left-hand side 
involves propagators in one direction only. The derivation in Eq. (3.32) shows 
that this is equal to the diagram in the right-hand side. 



in Figure 3.5 is calculated, which will show how the expressions for the T^ % ' can 
be derived. The derivation for channels (1) and (2) is completely analogous but 
involves the T^ ph ^ instead of T^ PP K The integral representation of this diagram is 
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Now the pole structure of the pp interaction (Eq. (3.20)) has to be inserted into 
Eq. (3.30) in such a way that the final result 



Eq. (3.30) 
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To arrive at an eigenvalue equation, we propose a discrete pole structure for 

■y(i)(2plh)n , v(2pUi)nu 
"p^hup'^h'^) (2plfc) . ' ^""J 

where in this equation 

y(2plh)n _ y{i){2plh)n ,~ „ A \ 

Substituting this form of i?« into Eq. (3.29) and multiplying by (E-eS' plh) ) while 
simultaneously taking the limit for E — > e^ 2pl/ ^ and assuming that the correlated 
i? has no poles in common with the free 2plh propagator G(°)> gives, again for 

W = (3), 

,,(3)(2plh)n = \p r (0)> ( (2plh) )r (3) ( Mplh)^ 

PlP 2 h! PiP2h' 1 ,pip 2 h 1 \ fc n / PiP2hi,p 3 Pih 1 \ t n I 

< x (l)(2plh)n x (2){2plh)n\ (3.35) 



P3,Pi 

X 



To go any further in this derivation, it is necessary to postulate a property of the 
pp vertex. When the integral in Eq. (3.14) is worked out and the limit procedure 
is applied, but now for E — > e a + ep, the result is a crucial property of the vertex, 
namely 

rS 75 (e a + e/3 ) - 0. (3.36) 

This property is of great importance to reduce the complexity of Eq. (3.35) together 
with Eq. (3.32). An analogous relation holds for the ph vertex T^ ph \ After some 
algebra the equations can be transformed into 



■y(3)(2plh)m 
P\Pih\ 



^ / jj(PP) n f?j(PP) n 



7 jxppyn j A 
. L/l PlP2 l M P3P4 
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( V (pp)n\* v (pp)n 
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1 tn t pi t p2 I I t n t p3 t P4 I 

This is a non-linear eigenvalue equation with the Faddeev energies em Pl ' 1 ' as eigen- 
values and the amplitudes of the Faddeev components as eigenvectors. Eq. (3.37) 
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allows us to define the following vectors 



p 1 p 2 h 1 ,nh' 1 1 1 y/2 

M . = ^>j=y^\ n (3.39) 

(t {3) ) , h , = (3-40) 

V / pip2hi,nft' 1 V2 v 7 

(# (3) ) _ ffci - (e^)+- e/tl ). (3.41) 

The definitions for (i) = (1), (2) are analogous but involve the eigenvalues and 
eigenvectors of the ph RPA. The second index of f/W, ifW, yC*) and both indices 
of £)W are a combination of an index labeling the RPA eigenvalues and the line 
propagating freely. This guarantees that both indices include the whole 2plh space. 
In these vectors, the formulation of the non-linear eigenvalue equation for any (i) 
is given by 

*(i)(2 P i/o = (Vw— -i — rWt + irWffWt 

x ^(j)(2pl/ l ) + ^(fc)(2plfc)^ _ (3 42) 

The index (i) can also be captured in a vector formulation by introducing 

( X (l)(2plh) \ 
AT(2)(2pifc) (3.43) 
X (3)(2plh) J 

and the matrix M in this space of three-component vectors 



M = 7 7. (3.44) 




These definitions together with the convention that the matrices U, H, T and 
D are the matrices with their respective components on the diagonal allow the 
reformulation of all three components of Eq. (3.42) into one matrix equation 

x = ( v^)-u rt + HH ") Mx - (3 - 45) 

After some algebra the problem then reduces to a linear non-Hermitian eigenvalue 
equation: 

luX = FX, F = A^B, (3.46) 
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with 



j _JT(i)jT(Dt -tfWtfMt \ 

A = | -HWHW I -HWH<M (3.47) 

.jy(3)_H-(3)t __ff(3)jy(3)t / y 

[/(i)y(i)t (/(i)r(i)t \ 
B = ! [/(2)T( 2 )t [/(2)r(2)t (3.4 8 ) 

[/(3) r (3)t [/(3) T (3)t q y 

£/(!)£,(!) £/(!)-! 

+ [ [7(2)^(2)^(2)-! o 

tf(3)£)(3)[/-(3)-l 

_ H (2) H m j _ ff (2) ff (2)t ] (3.49) 



3.5 Elimination of spurious solutions 

Through the definition of the three Faddeev components, the dimension of the 
matrices that have to be diagonalized has tripled. The number of physical solutions 
to the Dyson equation has stayed the same, so two thirds of these solutions have to 
be discarded. This can be seen from the way in which the final amplitudes X^ 2plh ^ 
are retrieved from the Faddeev components (Eq. (3.34)). Combinations of Faddeev 
components that sum to zero have a vanishing amplitude in the irreducible 2plh 
propagator and will not contribute to the selfenergy. 

The solutions to the three-component eigenvalue equation can be divided into two 
classes according to their eigenvalue with respect to the operator 



P = I ex , (3.50) 




where I ex is the matrix that exchanges the two particle indices 

7 P?P2hi,pip£hi = 5 PiP , 2 5 P2P' 1 5 h 1 h' 1 - (3.51) 
P commutes with F and the common set of eigenvectors of these operators are 



j exx (2pih) | and ,Y| 

,(2plh) Jex v (2plh) 



/ xi 2plh) 

je XX (2 P lh) | (3 52) 

,(2plh) . Texv (2plh) 



x \^p*-<<-) _ jex X W-"/ j y x \ 

with Xa 2plh ^ and xP plh ^ arbitrary vectors in the Faddeev component space. It 
can be shown that X_i and X\ are also eigenvectors of F. X_i gives rise to 
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anti-symmetric amplitudes to the irreducible 2plh propagator, while X\ results 
in symmetric amplitudes under the exchange of the two particle indices. Since 
we are dealing with electrons, the eigenvectors of interest are of the type 
They consist of the symmetric and anti-symmetric part of X^ 1 ^ and only the 



anti-symmetric part of X b 
Defining the bases 



(2plh) 



Now the unphysical solutions can be selected. 




and Y Sp = 




(3.53) 



it can immediately be seen that the vector space spanned by Y Sp concerns eigen- 
vectors that sum to zero amplitude and has to be projected out. Since we only 
need the physical solutions of the eigenvalue equation for the spectral amplitudes, 
it suffices to project both matrices A and B from the right-hand side with the 
physical and spurious solutions and take the matrix block (Y ph \F\Y ph ) to be 
diagonalized[Bar02, NBG99], 

The second row and second column of the projected matrix can be left out because 
of the exchange symmetry. The right-hand side projection of A results in 



A = 



A ( Y Ph Y s p ) 

J / I _ _ 2 #(l)ij(l)t (/ _ jes) 

71 V I~I ex - 2H^H^ (I - I ex ) 



(3.54) 



-i-hWhm 

I _jex + H (3) H (3)i (J_ je 
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71 



j _ jex _ 2/f (l)#(l)t (/ _ je*) 
I - jex _ 4Jj(3)_ H -(3)t 

V 3(/-/ ex ) -au^uw 

Q R 
P S 



-i-hWhm 

I - I ex + 2H^H^ 
2[/(3)t/(3)t 



(3.55) 
(3.56) 

(3.57) 

(3.58) 



where the transition between the last three lines in the equation are accomplished 
by using the anti-symmetry in the first two indices and by inserting the closure of 
the RPA equations Eq. (2.95) and Eq. (2.108). The projection of B goes along 
the same lines: 



B 



1 

71 



{-2U^U^Ef + SUWdWTJW- 1 ) (I - 



L 

M 



N 
K 



2f/3£/(3)t£/ 

(3.59) 
(3.60) 
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where the diagonal matrix was introduced as 

£ PiP2hi ) pipihi = + - e h) 5p 1 p> 1 5 P2 p> 2 Sh 1 h< 1 - (3.61) 
The matrix A can be inverted by using block-inversion[PTVF07] 

A- 1 = §) (3-62) 



(Q-RS^P) 1 -QRS- 1 
-S-ipQ S- 1 + S^PQRS- 1 



(3.63) 



The matrix block of interest is 

(y p ' l |F|F p ' 1 ) = QL + RM (3.64) 
= Q(L- RS^M) . (3.65) 

The calculation of Q can be avoided by solving the general eigenvalue equation 

p-i x (2pih) e (2 P ih) = r L _ RS -i M ^ x (2pih)_ (3.66) 

The left-hand side is given in terms of the Faddeev component matrices as 

Q-i = i (J - - 2^< 3 >t/( 3 )+ 

+ ?7 (3) f/ (3)t C/ (l)t-l C/ (l)-l( / _ /ex-, (3 6y) 

while the right-hand side is given by 

L-RS^M = U^D^UM- 1 -2U^U^E f 

+C/ (3) [7 (3)t f/ (l)t-l jD (l) C7 (l)-l (J _ jez) (3 68) 

This is the non-Hermitian eigenvalue problem that has to be solved and which 
results in the poles and amplitudes of the irreducible 2plh propagator in the present 
approximation. 



3.6 Coupling to the single-particle space 

After the calculation of the excitation energies e^ 2plh ^ and e^ 2hlp ^ together with 
their respective amplitudes, these excitations have to be coupled with the single- 
particle space. The coupling happens through the procedure introduced in Eq. (2.81), 
applied in the energy domain. Instead of building the selfenergy matrix explicitly 
from the Faddeev eigenvalues and eigenvectors and diagonalizing it afterwards, it 
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is preferable to reformulate the problem in the space of the single-particle states 
and the 2plh and 2hlp solutions. In this space the hamiltonian matrix becomes 

H = I V^piMt e (2 P ifc) o j (3.69) 

V"(2/ilp)t o 

The two blocks e^ 2plh ^ and e^ 2hlp ^ are diagonal in the space of the Faddeev solu- 
tions. Note that there is no explicit coupling between the 2plh and 2hlp states. 
Their mixing blocks are equal to zero and these spaces are only linked indirectly 
through the single-particle space. The two blocks V < - 2plh ^ and V i - 2hlp) take into 
account this coupling. These blocks are not simply the bare interactions. Due 
to the procedure introduced in Section 3.3, where the class of diagrams was re- 
stricted, we have neglected a class of diagrams. These are, however, not energy 
dependent and can be reintroduced a posteriori by replacing the bare interaction 
with the interaction leading to the description that sums all diagrams up to third 
order. As was derived by Schirmer et al.[SCW83], the interaction that gives the 
correct behavior is 

T/ (2plfc) _ T/ I I V Vah 1 ,h2h 3 Vh2h3,p 1 p2 

V a,p 1 p 2 h 1 ~ &hi,piP2 ' cy e , _|_ F , _ e _ e 

h 2 ,h 3 h2 hs Pl P2 

Vqp 3 ,pih 2 ^h 1 h 2 .,P2P3 _ Vap3 ,P2h 2 Vh 1 h,2,piP3 



P3,h 2 P3,h 2 



(3.70) 



An analogous expression can be found for the coupling between the single-particle 
states and 2hlp states. As the diagonal blocks are expressed in the basis of the 
physical states Y ph from Eq. (??), a sum has to be carried out over the 2plh 
indices to change to the basis of Faddeev solutions 

nr h) = e &y^P2H,. ^ 

Pi,P2,hi 

The diagonalization of the Hamiltonian matrix Eq. (3.69) has to be done iteratively. 
In order to resum all diagrams correctly up to third order, the Hartree-Fock selfen- 
ergy has to be calculated using the correlated density matrix. At every step of the 
iteration, a new density matrix is calculated from the eigenvectors of Eq. (3.69). 
Based on this density matrix, a new Hartree-Fock term is generated and inserted 
in the Hamiltonian matrix. This procedure has to be repeated until convergence. 
Note that the coupling blocks and diagonal Faddeev blocks stay the same during 
this iterative process. Therefore, they have to be calculated only once. 
This partial selfconsistency procedure is also used in the FRPA calculations for 
nuclei[BD02]. However the procedure here is somewhat more involved. In order to 



50 



3. Faddeev Random Phase Approximation 



include the correct description of the short range tensor correlations, a presumma- 
tion of diagrams is carried out. This is done by solving the Brueckner-Goldstone 
equations for the effective nucleon-nucleon interaction or G matrix[MS93a]. In 
this procedure, the Pauli operator excludes the states from the space to which 
the many body perturbation theory is applied. This prevents double countings of 
ladder-diagrams, while at the same time adding the intermediate states that are 
necessary to obtain correct results. This presummation of the Brueckner-Hartree- 
Fock propagator constitutes the main difference between the processes of calcu- 
lation for nuclei and molecules. If one wanted to apply the same routine as for 
molecules, one could use renormalized low-momentum interactions. These inter- 
actions have the short-range physics included and a prediagonalization of ladder 
diagrams is no longer needed. The low momentum interactions do however require 
the use of three-nucleon forces in the formalism. 



The following section presents the proof that when interactions at the TDA level 
of theory are used in the Faddeev procedure, the Faddeev TDA corresponds to 
ADC(3)[SCW83]. The single-particle selfenergy block and coupling matrices in the 
Hamiltonian Eq. (3.69) are identical to what is found in ADC(3). The question is 
whether the diagonal matrices reduce to the standard ADC(3) form after elimina- 
tion of the spurious solutions. The reduction can be done by setting the matrices 
that involve backward propagating amplitudes to zero, expressing the lack of such 
amplitudes in the TDA. The corresponding eigenvalue equation for the 2plh space 
is 



where the matrices f/W, T"W and £)W are now built up with the RPA quantities 
replaced by TDA quantities. If one does the same substitution of zero backward 
going amplitudes in the closure relations Eq. (2.95) and Eq. (2.108), the first two 
terms reduce to 
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(3.72) 




) 



PiP2h 1 ,p' 1 p' 2 h' 1 



X ( e Pi + e P2 — € hi) 
+ ( ^pip' 1 6Vp^ l! / l ' i p 2 , 



(3.73) 
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while the last two terms give 

( 2 [/(3)r( 3 )t + tf(3)£,(3)tf(3)-i) = b hlh ,(6 Plp[ 5 P2p ,-6 Plp ,5 P2p ,) 

V / piP2/ll,PiP 2 "l ^ 

X Opi + £ P2 

3 

+ 2$h 1 h[V PlP2tP > iP2 . (3.74) 

Together with the anti-symmetrization of both particle indices this yields for the 
eigenvalue problem 

e (2plh) x (?plh)n _ ^ (5^^ ((^p'^pj, - 5 pip'Jp2P , 1 ) ( e Pi + e P2 - € hi) 

p'^p^K 

~\~fih 1 h' 1 Vp 1 p2,p' 1 p' 2 + 3pip' 1 Vp' 2 h 1 ,h' 1 p 2 ~ 3pip' 2 Vp' 1 h 1 ,h' 1 p 2 
+fip2P' 2 Vp'i>il,h' 1 Pl ~ 3p2P' 1 Vp' 2 h 1 ,h' 1 p 1 ) ^p'^h^ ' (3-75) 

which is exactly the equation that is solved in the ADC(3)[SCW83]. In what 
follows, results obtained with ADC(3) will be labeled as FTDA, seeing that both 
methods are completely equivalent. 

There have been efforts to make a connection between Green's function theories and 
CC theories. The equivalence of the ring CC with doubles (CCD) and RPA[SHS08] 
shows that in some cases it is possible to find such a connection, although this 
correspondence holds only for the RPA with direct interaction and is based on 
the equality of the correlation energies. When the level of theory becomes higher, 
the equivalence relations become more and more intricate. It is possible to link an 
approximate CC method with doubles (CC2) to the ADC(2) matrix by replacing the 
CC2 ground-state amplitudes by those obtained in first-order perturbation theory 
and taking the Hermitian average[H05]. There is no obvious systematic mechanism 
for how to extend these correspondences to higher orders or more complex theories. 
Given the very involved diagrammatic content of FRPA, it should be no surprise 
that no equivalent CC theory has yet been found. However, hybrid methods using 
concepts of Green's function theory combined with coupled cluster theory have been 
devised[NS92] and have proven to give accurate results for small molecules[NS93]. 



Chapter 4 



Results 



In this chapter the method developed in the previous chapters is applied to some 
physical systems. First an application to atoms is discussed. Considering the 
spherically symmetric nature of closed-shell atoms, it is possible to exploit the 
symmetry to the maximal extent and to push the calculations towards very large 
basis sets. Theoretical extrapolation to the basis set limit makes precise comparison 
with experiment possible. 

The application to small molecules is another example where quantum mechanical 
correlation plays an important role in the description of the physical properties. The 
most common quantum chemistry methods like HF and density functional theory 
(DFT) fail to simultaneously describe correct ground-state energies and ionization 
energies, even at equilibrium geometry. First the equilibrium properties of these 
molecules are studied, after which the dissociation is studied in some detail. The 
physical implications of the separation of the atoms in a diatomic molecule are 
catastrophic for the RPA. Triplet instabilities prohibit the calculation of properties 
in the FRPA. 

The Hubbard molecule is used as a toy model to investigate the this further and 
to search for possible solutions. 

4.1 Accurate atomic calculations 

The first application of the FRPA to electronic structure was a calculation for 
neon[BVIMD07]. This calculation showed that the main ionization energies ob- 
tained with the FRPA are at least of the same level, and even somewhat better 
than those obtained with ADC(3). Two basis sets were employed: the augmen- 
ted correlation-consistent polarization valence triple-zeta (aug-cc-pVTZ) basis, a 
standard quantum-chemical Gaussian basis set and the HF+continuum, a numer- 
ical basis set based on Hartree-Fock and subsequent discretization of the con- 
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tinuum. Comparing results obtained in a finite basis set to experimental values 
may lead to faulty conclusions. Therefore, it is of importance to investigate the 
behavior of the results in the basis set limit, i.e. to systematically enlarge the basis 
set until the result converges. We considered a set of neutral atoms and ions with 
closed shells and performed calculations using Gaussian basis sets. The smallest 
systems (He, Be and Be 2+ ) were described in correlation-consistent polarization 
valence basis sets (cc-pVXZ) of increasing quality from double to quintuple zeta 
(X=2-5). For the larger systems, it is necessary to include core-valence functions 
to describe the correlation correctly. For these systems the correlation-consistent 
polarization core-valence (cc-pCVXZ) basis sets with increasing X were used. The 
basis sets for the charged species were obtained from the neutral-atom basis sets 
by rescaling according to 



where Z is the central charge of the nucleus and y is the total charge of the ion. 
All calculations were performed without frozen-core approximation, in a full active 
space. Since the FRPA calculations are done in a non-relativistic framework, the 
experimental values need to be corrected for relativistic effects. This has been 
studied in Ref. [Per04]. Since both ADC(3) and FRPA are of comparable level of 
theory, we assume that the relativistic corrections of ADC(3) are also applicable 
to the FRPA values. Both the FT DA and FRPA results were calculated with self- 
consistency on the Hartree-Fock diagram. This means that all results in this section 
are in fact equivalent to the FTDAc and FRPAc results defined in the next section. 

4.1.1 Convergence of the atomic calculations 

Table 4.1 shows total binding energies for Ne in the cc-pCVDZ basis and for Be in 
four increasingly large basis sets. The results obtained in FRPA are shown together 
with values calculated in the FTDA, CC with singles and doubles (CCSD) and FCI 
for comparison. For the Ne atom, all three methods are in agreement with the FCI 
result, showing only a small correction of less than 2 mH. The discrepancy from 
the exact result of FTDA is halved by going to FRPA. The Be atom represents 
a hard test case for electronic structure methods. Due to the near degeneracy 
of the 2s and 2p orbitals, its ground state is very badly described by a single 
Slater determinant. As a result the RPA will be driven to an instability in the 
S = 1 channel. For this channel alone, the RPA eigenvalues and eigenvectors are 
replaced with their TDA counterparts. All other channels are calculated at the 
RPA level of theory. The results show that the values calculated with FTDA and 
FRPA are very close. This means that in the case of the Be atom, it is not crucial 
to use RPA for the screening of the interactions. This should not be surprising as 
the RPA is not the optimal method for few-body methods since the violation of 




(4.1) 
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the Pauli-principle that is inherent to RPA plays a more important role in these 
systems. It is reassuring that the FRPA and FTDA results are close in this case, 
implying that one can safely use the FRPA even for smaller systems. Both Green's 
function methods are systematically off by 9 mH compared to CCSD and FCI. This 
means that the description of the excitations at the level of RPA is not enough to 
describe all correlations. More involved methods for the excitations could cure this 
behavior. 



E tot 


Ne 




Be 








cc-pCVDZ 


cc-pVDZ 


cc-pVTZ 


cc-pVQZ 


cc-pV5Z 


ADC(3)/FTDA 

FRPA 

CCSD 


-128.7191 
-128.7210 
-128.7211 


-14.6089 
-14.6084 
-14.6174 


-14.6154 
-14.6150 
-14.6236 


-14.6314 
-14.6310 
-14.6396 


-14.6375 
-14.6371 
-14.6457 


full CI 


-128.7225 


-14.6174 


-14.6238 


-14.6401 


-14.6463 



Table 4.1: Total binding energies (in Hartree) for Ne and Be obtained for cc-p(C)VXZ 
bases of different sizes. The results obtained with ADC(3) and FRPA (with 
self-consistency in the HF diagram) and with CCSD are compared to FCI 
calculations. 



The data in Table 4.2 involve the extrapolation of two results according to the 
formula 

E x =E co +AX-\ (4.2) 

where X is the number of zeta functions in the basis set and is the extrapolated 
basis set limit. This relation is known to yield good extrapolations for correlation 
energies[HRJ + 04]. Table 4.2 indicitates that the extrapolated result changes by 
less than 2 mH when going from extrapolation between X=T,Q to extrapolation 
between X=Q and 5. This leads us to decide that the current accuracy suffices 
to calculate the basis set limit and the deviation serves as an error estimate of the 
extrapolation. The convergence for the larger atom Mg is slower but the error does 
not exceed 10 mH. 

The convergence of ionization energies (lEs) and electron affinities (EAs) tends 
to be faster than the total binding energies because they represent differences of 
energies. Inaccuracies are comparable for the N±l and the N particle system and 
they tend to cancel out. This means that Eq. (4.2) still holds for lEs and EAs. This 
relation is not always right. When shake-up configurations are important, the basis 
set limit will not obey this relation. Also for smaller basis sets the behavior of the 
extrapolation may be prone to larger instabilities. An example of both possibilities is 
given in Table 4.3. The 3p orbit of Ar has a strong one-hole character and therefore 
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has smooth convergence behavior. The 3s orbit of Ar is a notable counterexample. 
Here the calculated values show an oscillatory behavior when increasing basis set 
size. This does not exclude smoother convergence when one goes beyond X=5, 
but this kind of calculations are beyond our current capabilities. In what follows 
we apply Eq. (4.2) and estimate the errors for the lEs to be of the order of 2 mH. 

4.1.2 Extrapolated results for simple atoms 

Table 4.4 presents the extrapolated ground-state energies for a series of simple 
atoms. The formula from Eq. (4.2) was applied using the data from the calculations 
for the X=Q and X=5 basis sets. For the lighter atoms (He and Be) the cc-pVXZ 
basis sets were used, while for the larger atoms the cc-pCVXZ basis sets were used. 
The experimental values have been corrected for relativistic effects along the lines 
of Ref. [Per04]. The data obtained at the CCSD level for He and Be 2+ can be 
regarded as FCI results. Still there is no guarantee that the extrapolation to the 
basis set limit has fully converged. This explains why these values can be below the 
empirical values. For these two species, we see that FRPA and FTDA lack about 
1 mH of the correlation energy, which corresponds to 2% of the total correlation 
energy. Be forms a special case as explained earlier. The bad description of the near 
degeneracy precludes the correct description with the Green's function methods. 
It should be noted that the error for the CCSD calculation is also largest for Be. 
This may lead to the conclusion that the basis set used in this calculation does 
not allow for the correct description of the correlation effects in this case. The 
FRPA describes at least 99% of the correlation for the larger atoms. Both FRPA 
and FTDA agree with the experimental results within the error estimate of the 
extrapolation procedure. It can generally be said that FRPA predicts the ground- 
state energies at least as well as FTDA (a rms = 9.5 mH). When the deviating 
results for Be are kept out of the calculation of <r rms , the FRPA (a rms — 3.4 mH) 
achieves even slightly better accuracy than CCSD (a rms = 4.2 mH). 
The extrapolated ionization energies are presented in Table 4.5. Again the differ- 
ence between the light (He and Be) atoms and larger atoms shows. The FTDA 
results are approximately 1 mH closer to the experimental results than the FRPA 
results. This is due to the bad few-body behavior of the RPA. The FRPA overes- 
timates the correlations, which results in ionization energies that are too low. The 
ionization energies of the larger atoms are systematically better reproduced by the 
FRPA than in the FTDA. This results in a a rms that is 3 mH lower for FRPA than 
for FTDA. The FTDA and FRPA are drastically better than the second-order res- 
ults, which is recovered when R is the free propagator. This indicates that one has 
to describe the selfenergy at least up to third order in perturbation theory correctly 
to account for the correct behavior of the ionization energies. The FTDA is the 
lowest level of theory that takes into account mixing of 2hlp states, which is re- 
quired to describe the higher ionization energies correctly. One such example is the 
3s orbital of Ar. In second order, there is one peak at an energy 36 mH away from 
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the experimental value. This peak has a spectral strength of 0.81. The second 
smaller peak carries a strength of 0.10 and lies at a higher energy. The inclusion 
of configuration mixing between 2hlp states shifts both energies downwards and 
redistributes the spectral strength of the main peak to 0.61, which corresponds 
better to the experimental value of 0.51. 
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It can be concluded that both FTDA and FRPA yield very similar results for the 
lightest atoms, while for the larger systems, going from TDA to RPA for the pp 
and ph channel leads to a systematic improvement of the ionization energies. The 
accuracy of the FRPA for small systems gives confidence that the method can be 
used even for few-body systems. The extrapolation to the basis set limit shows 
that the FRPA yields the same accuracy as in smaller basis sets and does not lead 
to overcorrelation. 



4.2 Molecular results 

After calculations for nuclei and atoms, the next logical step in the study of the 
FRPA are calculations for molecules. The molecules that were chosen all have a 
closed shell structure and their ground state is a spin singlet. Due to this singlet 
ground state, the results are independent of the third component of the spin. Spin 
angular momentum algebra can be applied to all matrices, which improves the com- 
putational scaling. The calculations were again done in Gaussian basis sets[BSE]. 
With the exception of H 2 0, all molecules studied here are linear molecules. This 
means that there is only one parameter fixing the external potential: the nuclear 
separation. The calculation is facilitated through the conservation of the rotational 
quantum number around the nuclear axis. The resulting reduction of matrix size 
yields a considerable speed-up of the calculations and a lower memory usage. The 
results are compared to CCSD(T) and to FCI or experiments where available. 

4.2.1 Ground-state properties 

Table 4.6 and Table 4.7 present the ground-state properties of a series of molecules 
which are ordered according to increasing correlation energy. These ground-state 
properties were obtained by performing a scan of different separation lengths (and 
angles in the case of H 2 0) around the equilibrium geometry for each molecule. 
A third-order polynomial was then fitted to these data points, after which the 
minimum was located. It is this minimum that is tabulated here. The calculations 
were performed in a cc-pVDZ basis set. The FCI results were calculated with 
the geometry from the FRPAc fit, which within the accuracy obtained here also 
corresponds to the CCSD(T) fit. The CCSD(T) ionization energies were derived by 
performing two calculations, one for the N particle system and one for the N — 1 
particle system, subsequently subtracting both results. The CCSD(T) results for 
H2 are equivalent to FCI. The FTDA and FRPA columns present the results of the 
procedures without performing the iterative calculation of the HF-diagram. This 
corresponds to the first iteration in the procedure for the calculation of FTDAc 
and FRPAc, which do include the Hartree-Fock diagram with the correlated density 
matrix. 
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The ground-state energies for Table 4.6 are very similar for FTDAc and FRPAc, 
differing by at most 4 mH. The deviations are such that the FRPAc is closer 
to the FCI results where available. The FRPAc results are generally close to the 
values obtained with CCSD(T), with a maximum deviation of 5 mH for H 2 0. The 
bond lengths obtained with FTDAc, FRPAc and CCSD(T) are of the same quality 
when compared to the experimental values. The Green's function methods are 
slightly better at describing the first ionization energies. The importance of the 
iterative calculation of the HF-diagram is clear from this table. The results are 
systematically driven closer to the reference values by adding this diagram, for 
both FTDA and FRPA. 

The deviations between FTDAc and FRPAc for the ground-state energies are of 
the order of 10 mH. The FRPAc ground-state energies are generally closer to 
the CCSD(T) results, which will be regarded as the reference here. Calculating 
the geometry of C2H2 forms a problem for the Green's function methods, while 
in the other cases they predict geometries close to the experiment. It should be 
emphasized that for these molecules the presence of the iterative HF-diagram plays 
an even greater role. For the calculation of C 2 H 2 , N 2 and C0 2 it is even the case 
that no minimum is found in the potential energy curve without this diagram. This 
does not mean that there are no results, only that the dissociation curve has the 
wrong shape. For the other two molecules, self consistency tends to adjust the 
results toward experiment. 

4.2.2 Ionization energies 

To make the comparison with previously published ADC(3) calculations[TS05], 
Table 4.8 has an extra column with these results. This column agrees well with the 
values obtained with the FTDAc. The small differences are already present in the 
reference HF calculations. The basis set used for this table was the augmented- 
cc-pVDZ (aug-cc-pVDZ) basis set. The FRPAc and FTDAc are of comparable 
accuracy. The lcr u level of N 2 is an exception, FRPA and FRPAc are off by about 
0.8 eV, indicating that in this case the RPA is driven close to instability for this 
channel, which affects the results negatively. 

4.2.3 Basis set convergence 

In Table 4.9 the results for different Gaussian basis sets are shown for the HF 
molecule. The differences in ionization energies between the basis sets with double 
zeta functions and those with triple zeta functions are of the order of 0.75 eV for 
the non-augmented and 0.25 eV for the augmented basis sets. The convergence 
behavior of the ground-state energies calculated with the FRPA and FRPAc is very 
similar to that of CCSD(T) and MP3. The convergence of the FRPAc is faster, 
again indicating that the self-consistency for the HF-diagram is desirable. The 
ground-state energy in the largest basis set aug-cc-pVTZ obtained with the FRPAc 
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is very close to the CCSD(T) result and outperforms the value from FRPA and 
MP3. 



4.2. Molecular results 



65 



Q. 
X 
LU 



Q 
in 
u 
u 



u 
< 

Q_ 

cr 



< 



< 



LO H 00 

b- ^ o 

<-l l> ,n 



to 



T-H 
CO CO 
1 — I b- 



i — I b— CO 

!B lO O 

>-! ^ co 

i— I O i— I 



o o 
b- b- 

b- 



i — I b— CO 
O lO O 
H N ,n 



O O) ID 

b- co i— i 

<-l b- ,n 



&3 



o 
co 



CO 

co 

oo 



io Ol O) 
co co oo 
oo co • 



CM b- CO 

co co b- 
oo co • 



co n ^ 
lo oo oo 
oo co _j 



i — i b- 

co co 
oo co 



B? i- 
^ cq 



I LO 
b- 

CN 



^ o co 

LO Ol CM 

=1 =1 CM 

O 

CO 

I 



LO CO ^ 

LO Ol CM 
CM CM -C 



o 

CO 

I 



co 

Ol 



O 
CO 

I 



CO b- 

LO Ol CM 

CM CM _J 



o 

CO 

I 



LO 
Ol 

cm co 



o 

CO 

I 



*2- 



b- CM 

T - 1 

°. co 

O 



co 

CM 

o 
o 



oo o 

CM CM 

CM Ol 

o o 
o 



oo co 

CM i — I LO 
N Ol ,n 



O 
O 



b- CO 
Ol LO 



co 

b- 



O 

o 



co o 

CM t-h b- 

=1 °5 lo 

O O 

o 



LO O 

b- O b- 

-1 l6 

O O 

o 



oo 

LO O 
Ol T-H 



t-h b- CM 

^ CO O Ol 

CM Ol i-H ^ 

CO CD 

b- 



CO CM CM 

co co o 

IN CI H 



CO 
b- 

I 



co i — i 

^ oo 



CO 

b- 

I 



O ^ CM LO 

CO O 
CM Ol i-H ^ 

CO O 
b- 



CM Ol 
CO o 

b- 

I 



O 
< 



b- 
o 



X 

CP 

m 



o 
x 



X 



O 

tn 



« » U 1) ,1, 
2 +* Q ^ S 

ir o u 
£ <-> 

"Oil 

-n cu iri 



i i 

bO - 

_£ bo < 
i_ c CL 

^ < o; 



y E q; 
o _ 



CU .:= 



N 
Q 
> 



TO 

-o 2 



■ ~n 
"a 



o 



O 

cu 
H 



HH 3 4-» 



o 
o 

CM 

o 



cu £= 
. o 



a. o 

s- ro 

ro -t-i 

in ™ 

5 ™ 



ro 



u a; 

O CL 
LL X 



■43 c o <u o 



bo CD 

= < 
5 H 



ro c o 

1 ™ a 

4-> < 

P Q 0- 

™ to 01 

>, U 

" u 



cu 



4-> O 

.2 - 

ui m 

c > 

o u "u 



o x 



co "D 

c O 

2 5 

' CD 



ro 



cn ro 

'o E 
.5 ' 



if) 

1- <L» 
J I 



1= 



U 

(U 3 (U 

S U h 



66 



4. Results 



o d 
3 o 

* 3 



n 
n 
n 



5 > 

73 ,5 

Tl <W 
> ft 

n 3 
g 3 

Q- dj 

p- 

0) — 



(0 



TO Tl 
O 7! 
C TJ 

§.> 



n 



era 
n>' 

S 

0) 

- ( 

n> 

n 



n 3 
dj_ qq_ 

o n> 
E_ 01 
dj -■ 

r+ 3 



w oq 

S 3 



QJ C_ 
r+ fD 



rt 2 



(/I J> 



rD n 

■n Q 



"3 <u 

fD 

i I 

ft o 

3 8 



r+ > 



rp 

Tl fD 

n 

fD c 

s g: 



SL < O i3 



<p_ fD 



5 



DJ " 
T3 fD 



Q_ zr 

fD 



QJ 

1 9 

§ < 
Q_ O 



E. rt 

» § 

fD =J 

Q_ " 

a. S 

DJ 3" 

r+ — 

DJ fD 



o 



'30 

oc 



oo u u u 



to a> co 
Or to CO 



oo 

00 



OS CO 
bO tO J> 



oo 
oo 



-J OS 

GO tO 



co 
-J 



o 

CO 

to 



o 
h co 

h- 1 to 

~j 

01 tO 



CO ^ ^ Cn 
io s ^ c 



o 

CO 

to 
-J 



. o 

cn co 

oo oo 



o 
o 



o 

to 

M 



a 



w U b 
GUM 

co o 



4^ • 

. I— 1 I— 1 

to CO o 

co co o 



tt^ to ^ 
co oo 



M 1 — 1 CO 

" ^ b 

Oi Ol 

^ Ol Cn 



o to 

oo 



^3 



to 

J> 

S -1 to co 
co oo co 
Cn Cn H 



I 

-a 

I — 1 M 1 — 1 -<t 

b to f-> 
to oo co o 
a> co oo to 



I 

i — 1 i — 1 i — 1 -<r 

t -1 b to b 
h oo co co 
^ o oo co 



I 

-a 

I — 1 M I— 1 ~-^J 

i -1 b to f-> 
o oo co i — 1 

oo H (O M 



. o to 

^ m o 
co co co 



to 

J> 

P to co 
s oo oi 
Cn 4^ Cn 



to 

co co 
o co 



to 

J> 

P to co 
co oo Oj 
Cn 00 



to 

isj CO 
co oo 

Cn O 



to 

~4 



4.2. 



Molecular results 



67 



4-J 

Q_ 


LO 
O 


o 


T— 1 

o 


T— 1 

Ol 


CM 


O 

ID 


co 

Ol 


CO 

1"-- 


CM 
ID 


h- 


i— f 

LO 


Ex 


id 
i— i 


o 

CM 


<* 

T— 1 


id 
i— i 


oi 

i— 1 


LO 

I— 1 


id 

T— 1 


cd 
i— i 


CM 
i — 1 


«^ 
i — i 


co 

T— 1 



u 



u 

< 



u 
< 
Q 
h- 



< 
Q 
I- 



LO 




Ol 




Ol 


CO 


^- 


o 


r- 


oo 


oo 


oo 


CM 


iq 


co 


LO 


I— 1 


I— 1 


Ol 


iq 


Ol 


i— i 


id 


O 


oo 


id 


oi 


LO 


h- 




CM 


<^ 


oi 


i — i 


CM 


i — i 


i— i 


1 — 1 


T— 1 


T— 1 


i — i 


i— 1 


I— 1 


T— 1 


LO 


00 




LO 


ID 


ID 




Ol 


CM 


1 — 1 


LD 


o 


O 


oo 


Ol 


<* 




h- 


CM 


LD 


Ol 


O 


id 


CD 




id 


oi 


LO 


h- 


CO 


CM 




O) 


i — i 


CM 


T— 1 


i— i 


1 — 1 


I— 1 


i— 1 


t— 1 


i — l 


1 — 1 


T— 1 



oo 


ID 




co 


Ol 


CM 


LO 


ID 


ID 


LO 


i— 1 




00 


Ol 


Ol 


1 — 1 


r- 


co 


o 


00 


i — 1 


CM 


lO 


o 


oo 


LD 


o 


LO 


ID 


Ol 


CM 


LO 


Ol 


1 — 1 


CM 


T— 1 


i— 1 


CM 


I— 1 


I— 1 


T— 1 


i — 1 


i — 1 


T— 1 



ID 
<* 


oo 
oo 


co 

00 


oo 

Ol 


i — l 
i— 1 


LO 

ID 


CM 

co 


Ol 
Ol 


oo 
oo 


1 — 1 
1 — 1 


Ol 
I— 1 


ID 
i — l 


o 

CM 


oo 
i — i 


LD 
i— l 


o 

CM 


LO 

T— 1 


ID 

T— 1 


CO 

T— 1 


CM 
I — 1 


LO 
1 — 1 


oi 

T— 1 



CM 
CM 


T— 1 


co 


CM 

o 


LO 

o 


i— l 


O 

CM 


LO 

OO 


o 

00 


ID 
O 


LO 
i— 1 


id 

1 — 1 


o 

CM 


1 — 1 


i— i 


CD 
CM 


ID 
i— l 


h- 
i— 1 


oi 
I— 1 


CM 
i — 1 


LO 
i — 1 


oi 

T— 1 



i — i 


00 
Ol 


o 

1 — 1 


**■ 


Ol 
Ol 


LO 
CM 


oo 

h- 


LO 

CM 


ID 
00 


oo 

Ol 


ID 

LO 


i — i 


o 

CM 


LO 

I— 1 


t— i 


I— i 
CM 




ID 
i— l 


T— i 

CM 


oo 

I — 1 


LO 
1 — 1 


oi 
I— 1 



t= b 
^ — i oo 



b t= b 

Ifl H 



b t= b 

OO i— I CM 



i — i oo 1 — i 



fa 
X 



O 

O 



o 

X 



ID CM 
CM ID 



i — I CO 

oo oo 



oo oo 

CM 



LO 00 
CM h- 



O O 

oo 



o o 

00 



00 00 
CM ID 



CO 
CM ID 



i — I ^ 
00 ID 



00 ID 



CM 



ID 
CM 



> A 
CD 

- — - H 
l< I 



68 



4. Results 



Method 


Level 


cc-pVDZ 


aug-cc-pVDZ 


cc-pVTZ 


aug-cc-pVTZ 


Expt. 


FRPA 
















Eo 


-100.172 


-100.106 


-100.335 


-100.305 


- 




Itt 


15.46 


16.05 


16.19 


16.33 


16.11 




3(7 


19.56 


20.03 


20.05 


20.22 


20.00 


FRPAc 
















Eq 


-100.228 


-100.261 


-100.346 


-100.357 






Itt 


15.54 


15.35 


16.16 


16.41 


16.11 




3(7 


19.54 


20.24 


20.00 


20.27 


20.00 


CCSD(T) 
















Eo 


-100.228 


-100.264 


-100.338 


-100.350 






Itt 


15.44 


16.06 


15.96 


16.16 


16.11 


MP3 
















Eo 


-100.224 


-100.256 


-100.330 


-100.340 






In 


15.42 


15.99 


15.88 


16.04 


16.11 



Table 4.9: Ground-state energies in Hartree and vertical ionization energies in electron- 
volt for HF, calculated in different basis sets. The geometry was taken at 
the experimental value of 0.917 Angstrom. Experimental values are from 
Ref. [TS05, KKA+81], CCSD(T) and MP3 calculations were performed at 
the same geometry. 



4.2.4 Example of the pole structure 

Figure 4.1 presents a plot of the spectral strength of the Green's function for 
the negative energy region of HF, calculated in the cc-pVDZ basis set. Only 
peaks with a spectral strength higher than 0.05 are retained. This plot shows the 
fragmentation of the deeper lying peak corresponding with the core orbitals of F. 
There are two peaks close to the Fermi-level that are almost unaltered compared to 
Hartree-Fock. Their spectral strength is around 95% in both FRPAc and FTDAc. 
The main difference with Hartree-Fock occurs for the peak around —40 cV. This 
peak is split up into two separate peaks with reduced strength. For FRPA this is 
a peak of 0.60 and one of 0.30, while FTDA gives a peak of 0.50 and 0.40. There 
are also some satellites at lower energies. 

The calculations at equilibrium geometry indicate that FRPAc yields results that 
are of similar quality as those obtained with FTDAc, just as for atoms. While 
RPA was developed for extended systems, it seems that FRPA is can be applied 
to both small and more extended molecules. The self-consistent treatment of the 
HF-diagram is mandatory to obtain good results. 
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Figure 4.1: The pole structure of the Green's function for HF in Hartree-Fock, FTDAc 
and FRPAc. The energy axis is in eV. 
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A point that deserves some attention is the conservation of particle number as 
predicted by the Baym-Kadanoff theorem. Since we are using only partial self 
consistency, it is to be expected that the particle number will not be exactly equal 
to the physical value. This can be seen from Figure 4.2, where the deviation 
from the real particle number is plotted in function of the separation distance of 
the nuclei. The deviation is of the order of 10~ 3 and is not varying much in the 
neighborhood of the equilibrium. The deviations are more strongly oscillating in 
the dissociation limit where the behavior of the RPA is getting worse. Surprisingly, 
there is almost no difference between FRPA and FRPAc in this respect and hence 
only the FRPAc results are shown in the figure. 



0.004 



0.002 



-0.002 



-0.004 




0.85 0.9 0.95 1 1.05 1.1 1.15 1.2 1.25 1.3 

r 

Figure 4.2: The deviation of the particle number for HF in the cc-pVDZ basis set. The 
deviation AiV is given relative to the physical particle number and the bond 
length r is in Angstrom. 



4.3 Dissociation limit 

For chemical purposes, it is interesting to look at the dissociation behavior of a 
method. In the disscociation limit, the constituents of the molecule are placed far 
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apart. A size-extensive method yields the same result in this limit as when the 
two constituents are calculated separately. Figure 4.3 shows the dissociation curve 
of HF calculated with FRPAc and FTDAc in the cc-pVDZ basis set. Around the 
experimental equilibrium geometry, the FRPA and FTDA results are very similar as 
seen from Table 4.6. At larger bond lengths, however, the FRPA energy starts to 
diverge from the FTDA result. At a bond length of 1.287 A, the FRPA result even 
stops to give a result. This is not due to a numerical instability in the calculation 
of the e ( 2 P 1 ' 1 ) or e ( - 2hlp ^ but to an instability in the ph RPA. The closer to this point 
of instability of the RPA, the less the RPA describes what physically happens. As 
a result the FRPA, that heavily depends on the quality of the intermediary states, 
starts to diverge and predicts wrong ground-state energies. FTDA continues to 
give qualitatively correct results. This is a trend that occurs for every molecule in 
the test set. Moving away from equilibrium, there is always an RPA instability at 
some point that prohibits the calculation of FRPA values. 
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Figure 4.3: The dissociation curve for HF in the cc-pVDZ basis set. The total energy 
Eq is in Hartree and the bond length r is in Angstrom. 



This behavior can be understood from calculations for H2 in a minimal basis set 
and is known as the "triplet instability" [Tho61, vP67]. It is not so much a problem 
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of RPA as of the HF wave function. In the dissociation limit, the HF ground 
state becomes unstable with respect to ph excitations. The Hartree-Fock stability 
matrix is nothing more than the matrix in the left-hand side of Eq. (2.89). When 
one of the eigenvalues of this matrix becomes zero, the RPA eigenvalues reach 
a bifurcation and become complex. The inadequacy of the HF ground state can 
be easily seen from an analysis of H2 in a minimal basis set. The spatial wave 
functions are Is functions centered on the H-atoms A and B. These can be put 
in a bonding and anti-bonding combination for the molecular orbitals. This is also 
the definition of the particle and hole state of the Hartree-Fock solutions 



\h) = \b) = ±{\A) + \B)) 



\P) = \a) 



V2 



(\A)-\B)). 



(4.3) 
(4.4) 



These states are normalized to one at great separation, which is the case of interest 
to us. The restricted Hartree-Fock ground state is the spatially symmetric state 
with positive parity and spin S = 



1*0} = 



1 



M s =0 



lit)) 

|0>. 



(4.5) 
(4.6) 



Particle-hole excitations can only be formed by removing a particle from a bonding 
state and replacing it with a particle in an anti-bonding state. The excitations can 
be split up in a spin triplet and a spin singlet. The energy of the triplet state 



|d>) 



<4 ® a\ 



s=i 



M s 



|0> 



(4.7) 



is found to be 



<$|H|$) = (ab\H\ab) - (ab\H\ba). (4.8) 

When the A and B H-atoms are separated by a large distance, the overlap between 
the wave functions becomes neglegible and the energy of the triplet state reduces 
to 



($\H\$) = (AB\H\AB) - (AB\H\BA). (4.9) 

This is the energy one would expect for the dissociation state where each H-atom 
receives one electron. It is the physically correct solution which restricted HF fails 
to describe. So the energy of the triplet state will always be lower than that of the 
restricted HF singlet ground state. 
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As a result of the instability of the HF ground state to a ph excitation, the RPA 
matrix no longer is positive semi-definite. As the separation of the two atoms grows, 
the lowest S = 1 RPA eigenvalue approaches zero and then becomes complex. This 
behavior is plotted in Figure 4.5. At a separation of 1.2 A, the excitation energy 
crosses zero and becomes complex. As a result, there are no FRPAc results beyond 
this separation distance in Figure 4.4. The TDA does not have this problem. 
Because of the unconnected positive and negative energy diagonalization problem, 
a TDA eigenvalue can always be found. However, Figure 4.5 shows that the 
TDA eigenvalue in the same channel also becomes negative. This means that the 
Hartree-Fock state should not be used as a reference. It is possible to calculate 
FTDA values beyond this point even though the HF ground state becomes an 
invalid starting point for perturbation theory. 




Figure 4.4: The dissociation curve for H2 in the STO-6G basis set. The total energy 
Eo is in Hartree and the bond length r is in Angstrom. 
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Figure 4.5: The lowest S = 1 excitation energy as a function of the nuclear separation 
distance r for H2 in the STO-6G basis set. The excitation energy e^ ph ^ is in 
Hartree and the bond length r is in Angstrom. 



4.4 The Hubbard molecule 

In condensed matter physics, many efforts have been made to account for the 
behavior of spin systems. They form a model for the great variety of types 
of magnetic arrangements that occur in experiments with real solids. Since its 
introduction[Hub63], the Hubbard model in particular has been the subject of ex- 
tensive research. Although it was originally postulated for the study of correlations 
of d-electrons in transition metals, it is often used for the description of effects 
such as itinerant magnetism and metal-insulator transitions. At the same time, the 
repulsive Hubbard model for two sites ("Hubbard molecule") can be used to study 
the competition between localization and derealization of electrons in a schematic 
model. It is this case of strong correlations where the RPA fails. This tool has 
been used to cure some of the deficiencies of the GW approximation[RGR09] and 
to study an approach that resums pp and ph interactions[RBR12] with a screened 
interaction. The physical interpretation, together with the access to analytical res- 
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ults, makes the Hubbard molecule the perfect background to study the dissociation 
behavior of the FRPA. 

The Hubbard molecule is a special case of the general Hubbard Hamiltonian for 
two sites. The Hamiltonian in the site basis is given by 



where i and j run over the two sites 1 and 2 and Ui is the third component of the 
spin with possibilities j" and \.. Note that in the repulsive Hubbard model both U 
and t are positive. It is advantageous to go over to states with a definite parity 
and third component of the spin by applying the rotation 

a la = (a+„ +a_ CT ) (4.11) 

a 2a = (a +CT - a_ CT ) . (4.12) 

The matrix elements of the interaction conserve parity and spin projection, which 
gives the interaction a block-diagonal structure. The matrix of the single-particle 
part of the Hamiltonian is even diagonal in this basis. We are interested in the 
H2 molecule-like behavior of the Hamiltonian, so there are two electrons in the 
system. This corresponds to the Hubbard model at half filling, two electrons for 
four states. In this sector, the exact wave function is given by 

l*S> = 1 == _(ETl+t + +) 

V(4t- Vl6t 2 + U 2 ) +U 2 

+ (it - y/m? + Tp} |-t-l>) , (4-13) 

which has an exact N = 2 ground-state energy of 

K = * + £-^±* (4,4, 

For the exact Green's function, it is important to have the difference of the ground- 
state energy with the energy of the excited states. The N + 1 and N — 1 systems 
are the 3 and 1 electron systems. The two-site Hubbard model with only one 
electron is a trivial problem. The ground state is the positive parity state with spin 
projection f or |. This degeneracy has to be lifted by choosing one of the two 
projections. The excited states are listed in Table 4.10. 

The three-electron system forms the dual problem: one hole has to be located in 
one of four states. In contrast to the N = 1 system, the ground state here is 
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i 


*i 




1 


l+t) 


e - t 


2 


l+i> 


e -t 


3 


l-t) 


ea+t 


4 







Table 4.10: Ground and excited states of the Hubbard molecule in the 1 electron sector 



formed by a negative parity state. The hole state that has to be accommodated 
gets a minus sign for its parity so that it is a positive parity state. The excited 
states for N — 3 are given in Table 4.11. 



i 


^3 

I 






E? 


1 


l+t + 4-- 


t) 


3e 


-t + U 


2 


l+t + 4-- 


1) 


3e 


-t + u 


3 


I-1--4.+ 


t) 


3e 


+ t + u 


4 


I-T-I + 


1) 


3e 


+ t + u 



Table 4.11: Ground and excited states of the Hubbard molecule in the 3 electron sector 



After close inspection, it can be seen that the Green's function is diagonal in both 
parity and spin. The four only components of the Green's function are defined by 



1 " " L. E -{El-El) + l r 1 + ^E-{Et-El)~ lV ^ ] 



The numerators are the squares of the coefficients in the definition of the ground- 
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state wave function in Eq. (4.13). The resulting Green's functions are 

(4t- VlGt 2 + U 2 Y 



G+ a — 



G- 



+ 7 — ; I ( 416 ) 

1 




(4t-Vl6i 2 + t/ 2 ) 2 + (7 2 \g- (e -t+g + Vm ^ U * )+ir ] 
+ (^^F^f , . (4 , 7) 



IT] 



The resulting pole energies are plotted in Figure 4.6. In the limit of infinite interac- 
tion strength y — >• oo, the forward propagating pole exhibits linear behavior, while 
the backward propagating pole converges to U . Figure 4.7 shows the correspond- 
ing spectral strengths. At small interaction, only the backward propagating pole 
carries any strength. In the limit of higher interaction, both poles converge to a 
spectral strength of 0.5. This means that in this limit, they are equally important 
and the system can only be described by a fragmented single-particle propagator. 

To build the RPA matrices, it is necessary to dispose of the Hartree-Fock solution 
for the N = 2 system. The hole states are given by the positive parity states, while 
the negative parity states are found to be the particle states. The HF energies are 



e h = e o-t+^ (4.18) 



for the hole states and 



eo+t+^ (4.19) 



for the particle states. With the HF reference state, the ph RPA matrix (Eq. (2.89)) 
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Figure 4.6: The pole energies for the exact G+ CT are shown as a function of the reduced 
interaction strength y. The offset eo is chosen to be zero. The red line 
represents the forward propagating pole, while the blue line represents the 
backward propagating pole. 
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Z 




Figure 4.7: The spectral strength of the poles from Figure 4.6 are plotted as a function 
of the reduced interaction strength ~ for G+ CT . The red line represents 
the forward propagating pole, while the blue line represents the backward 
propagating pole. Both strenghts converge to 0.5 in the limit of infinite 
interaction. The pole strengths of are the same, but with the forward 
and backward propagating pole switched. 
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that needs to be diagonalized is 



A 

-B* 



( 2t 



u 

2 





u 

2 



B 

-A* 



2*-f 



u 

2 





2t 

u 

2 



U 
2 

2t 



u 

2 



-2t+% 




2 





-2t + 





U 2 



2 U 





u 

2 



-2i 

_ u 

2 



J7 
2 





u 

2 

-2t / 



The eigenvalues of this matrix are 



= + 2tU 



Aph)< 



^ h)> = V^ 2 - 2tU 
= -\/4t 2 - 2tU, 



iph)< 



(4.20) 

(4.21) 
(4.22) 
(4.23) 
(4.24) 



where the solutions are ordered according to their total spin in a singlet and a 
triplet. It is clear that for the triplet, an imaginary eigenvalue is possible. When 
the interaction y is greater than two, the eigenvalue becomes complex. This does 
not manifest itself in the case of the TDA eigenvalues 





(4.25) 


£ ° _ M 2 


(4.26) 


e<* h)> =2t-^ 


(4.27) 




(4.28) 



but at an interaction j of four, the excitation energy for the forward propagating 
pole becomes negative. This means that the excited state is lower in energy than 
the ground state of the system. This situation is illustrated in Figure 4.8. For 
reference the exact excitations to the N particle system are given in black. It is 
clear that the lowest two exact excitations are reproduced in the TDA and RPA. 
The spin singlet has no instability, while the RPA spin triplet becomes complex 
at an interaction j = 2 and the TDA triplet becomes negative at an interaction 
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j- = 4. The highest excitation has no counterpart in RPA or TDA because these 
two models do not allow the correct parity. The pp RPA channel does not exhibit 
this kind of critical behavior. Over the whole parameter range, the eigenvalues stay 
well behaved. 

The pp and ph RPA can now be used as building blocks for the FRPA. The FRPA 
matrices get too involved to be calculated analytically. The results of the simulation 
can be seen in Figure 4.9. The ground-state energy of the Hubbard molecule is 
plotted as a function of the reduced interaction strength. Just as with the real 
H 2 molecule, the results of the FRPA follow the exact values very closely until the 
breakdown point at j- = 2. So unlike the HF molecule, the results are trustworthy 
even for this system with a very strong few-body character. After the instability 
occurs, the calculation becomes impossible. The exact results coincide with the 
FTDA results. Even the poles and strengths are reproduced exactly. So for this 
problem, although the excitation energies are not exact in the TDA (see Figure 4.8), 
the FTDA Green's function is the exact result. 

One solution for this problem is to replace the unstable channel with its TDA coun- 
terpart, i.e. replacing one single energy and amplitude with the values calculated 
from a separate TDA calculation. This allows for values to be calculated beyond 
the instability, but the results obtained there are very close to FTDA. So no real 
advantage is gained. Another idea is to see if added self-consistency can cure this 
problem, as was done for example in Ref. [BD02]. This approach, however, does 
not eliminate the possibility of complex RPA eigenvalues. To accomplish a higher 
level of consistency, we use the exact Green's function to construct the II' ) in 
Eq. (2.84). Seeing the small dimensions of this problem, the diagonal nature of 
the exact Green's function and the division of the matrices into blocks with defin- 
ite spin, it is possible to directly invert Eq. (2.84) to obtain the RPA polarization 
propagator II. The resulting poles are shown in Figure 4.10 and Figure 4.11. A 
striking difference is the occurence of two times two branches per spin. This is 
due to the fragmentation of the exact single-particle Green's function, this doubles 
the matrix dimension of the RPA problem. Of course it is still impossible to reach 
the excitation with different parity, so there are two excitation energies too many. 
We consider the higher two excitations to be spurious and discard them in further 
calculations. It should be noted that the S — excitation becomes exact up to 
numerical errors. The 5 = 1 solution differs from the exact result, but it is striking 
that no instability appears at any interaction strength due to the use of a fragmen- 
ted single-particle propagator. Therefore, it is possible to do calculations beyond 
Y = 2 with the FRPA using these eigenvalues and eigenvectors. 

Figure 4.12 shows the results of this calculation. At low interactions, just as with 
the normal FRPA, this calculation based on the exact single-particle propagator 
follows the exact result very closely. Even up to a reduced interaction of four, 
the general trend of the exact result is reproduced. Afterwards the values become 
detached from the exact curve very quickly. At y = 10, the system is not even 
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Figure 4.8: The ph excitation energies for the Hubbard molecule are plotted. The red 
curves correspond with the S = 1 solution, while the blue curves are the 
5 = solutions. Dashed lines are TDA, full lines are RPA. The black lines 
are the exact excitation energies. It should be noted that the upper exact 
solution has no counterpart in TDA nor RPA. This is because the parity of 
this excitation is negative, while this impossible to reach with a TDA or an 
RPA ph excitation which removes a hole with positive parity and adds a 
particle with negative parity to the HF ground state. 



4.4. The Hubbard molecule 
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Figure 4.9: Ground-state energies of the Hubbard molecule versus reduced interaction 
strength. The red squares are the FRPA results with self-consistency on the 
level of the HF diagram; the black line is the exact result which coincides 
with the FT DA result. At ^ = 2 the RPA-instability sets in. 
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Figure 4.10: The ph excitation energies for the Hubbard molecule calculated with the 
exact propagator are plotted for the S = channel. The squares are the 
values calculated with the exact propagator as a reference, while the solid 
red line are the normal RPA values. The black lines give the exact S — 
excitation. The highest exact excitation is the positive parity excitation 
which is not produced in either the normal RPA or the RPA based on the 
exact propagator. 
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Figure 4.11: The ph excitation energies for the Hubbard molecule calculated with the 
exact propagator are plotted for the S = 1 channel. The squares are the 
values calculated with the exact propagator as a reference, while the solid 
red line are the normal RPA values. The black line gives the exact S = 1 
excitation. 
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Figure 4.12: Ground-state energies of the Hubbard molecule versus reduced interaction 
strength. The red squares are the FRPA results based on the ph RPA with 
the exact propagator with self-consistency on the level of the HF diagram, 
the black line shows the exact result. There is no instability in this case. 



bound anymore. It is unlikely that adding further self-consistency would positively 
influence the results, as the RPA is already calculated with the best propagator 
possible. This could mean that the RPA misses a crucial diagram that cannot be 
cured by self-consistency, but only by making the transition to a ph interaction 
that is exact beyond first order in perturbation theory. 

Schuck and coworkers[SE73, DS90] have developed the so-called SCRPA that goes 
beyond the first-order RPA. The SCRPA is a theory that can be applied to both the 
pp and ph interactions and differs from the standard RPA mainly in the use of a cor- 
related ground state for the calculation of its expectation values. This means that 
the first-order approach in perturbation theory of RPA is upgraded to a nonperturb- 
ative variational theory, although not of the Rayleigh-Ritz type. The SCRPA has 
been applied to some test models and also to the Hubbard model[SS99, JSDB05]. 



4.4. The Hubbard molecule 
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In Ref.[JSDB05] it is found that the SCRPA solves the two-body problem exactly, 
while conserving the matrix size and characteristic structure of the RPA eigenvalue 
problem. Due to the flexible nature of the FRPA, the SCRPA holds the promise of 
making the FRPA exact for two-body systems and allowing the calculations to go 
beyond the point where the classical RPA experiences unstable behavior. However, 
the SCRPA has not been applied to realistic systems yet because of the appearance 
of the single-particle density matrix in the equations. The RPA amplitudes give ac- 
cess to only part of the density matrix, while all matrix elements are needed to build 
the SCRPA matrices. This problem could be solved by incorporating the SCRPA 
into the FRPA framework. As it is a single-particle Green's function method, the 
FRPA gives easy access to the single-particle density matrix. A procedure could 
be developed in which the FRPA density matrix is used to update the SCRPA 
matrices, after which the SCRPA amplitudes and eigenvalues are used to update 
the FRPA matrix, and to repeat this until convergence. Taking a suitable starting 
point for the single-particle density matrix to be used in the first RPA iteration 
could circumvent the occurence of instabilities. 



Chapter 5 

Conclusion 



In this work we have introduced a Green's function technique for the calculation of 
ground-state energies and ionization energies in quantum many-body systems. The 
elementary building blocks of the theory are the RPA (random phase approximation) 
excitations. These are used to construct an approximation for the selfenergy in a 
specific manner, called the FRPA (Faddeev RPA). The FRPA holds the promise 
to be a Green's function method with a wide applicability, from finite systems like 
atoms and molecules to extended systems like the uniform electron gas and nuclear 
matter. The fully self-consistent FRPA is conserving in the Baym-Kadanoff sense 
and consequently obeys important conservation laws. 

In Chapter 3 we have presented the FRPA technique as a Green's function method 
that is able to couple pp (particle-particle) and ph (particle-hole) excitations of 
the RPA type. The procedure involves the extension of the space in which the 
selfenergy is expanded from the single-particle space, to the single-particle space 
together with the 2plh (two-particle-one-hole) and 2hlp (two-hole-one-particle) 
spaces. The consecutive inclusion of ph and pp excitations is governed by the 
Faddeev equations. The eigenvalues and amplitude from the RPA problems form 
the building blocks for the Faddeev components. The Faddeev technique is needed 
when including excitations at the RPA level. In this respect, the current method 
differs from the ADC(3) (algebraic diagrammatic construction method of third 
order) which only involves excitations of the TDA (Tamm-Dancoff approximation) 
level to be exchanged. Both methods are closely related: the ADC(3) can be 
rederived from the FRPA by the exclusion of backward going diagrams in the pp 
and ph interaction vertices. Diagrammatically, this close relation was already clear, 
but now it has been shown rigorously. Replacing the TDA intermediary excitations 
with their RPA counterparts implies that the FRPA behaves better in the limit of 
extended systems. It also means that the matrix equations become non-Hermitian, 
which can lead to the occurence of complex eigenvalues. This should always be 
kept in mind when performing an FRPA calculation. At first sight, the matrix 
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dimension of the FRPA seems to triple compared to ADC(3), but after carefully 
studying the nature of the solutions, we can discard two thirds as unphysical. The 
separation of physical and spurious solutions can be carried out by a straightforward 
projection, which reduces the matrix problem to its orginal size. 

First we have performed calculations of the electronic structure for a series of 
atoms. These calculations showed that even for very small systems, where the 
RPA violations of the Pauli principle are expected to be important, the FRPA 
obtained satisfactory results. In general, FRPA and FTDA yield very similar results 
for the lightest systems, while the additional diagrams in RPA lead to a small but 
systematic improvement as the atomic number increases. Both Green's function 
methods agree well with CCSD (coupled cluster with singles and doubles) in the 
basis set limit and deviations from experimental results are within the estimated 
extrapolation error. The near-degeneracy of the Be atom seems harder to describe, 
although the errors are systematic, appearing for all methods, and seem to be 
mostly due to deficiencies in the basis set. The conclusions for the ionization 
energies follow the same trend. For the smaller two-electron atoms, the FRPA 
does not improve on the FTDA results. For all other atoms, the FRPA pushes the 
calculated values in the direction of the experiment. The Be atom also forms a 
problematic case for the ionization energies. 

The next logical step is to calculate properties of small molecules. The diatomic 
molecules of interest have the nuclear separation as an extra parameter. Around 
equilibrium geometry, we find results that are in line with the atomic calculations. 
For smaller molecules, FRPA and FTDA are very close to each other, while for 
larger molecules the differences are more pronounced. In any case, the scan of 
nuclear separations clearly shows that self-consistency on the Hartree-Fock diagram 
should always be performed. The ionization energies are of comparable quality for 
FTDA and FRPA. In these calculations, there is also a problematic case: the 
third ionization energy of the N2 molecule produces results that deviate from the 
experiment. An attempt was made to push these results to the basis set limit, but 
this is for now still out of reach of our computer resources. 

The molecular calculations in the dissociation limit show the "Achilles' heel" of the 
FRPA. Being highly dependent on the quality of the RPA solution, the applicability 
of the FRPA only stretches so far as the applicability of the RPA. It is known 
that the ph RPA exhibits triplet instabilities when the bond in diatomic molecules 
is stretched. For the smaller systems, we do not find that the FRPA results are 
adversely influenced by the imminent RPA instability. However, the results for the 
HF molecule show that, for molecules with higher correlation energies, the results 
of the FRPA start to deviate in the vicinity of the instability. This means that the 
FRPA values can only be trusted around the equilibrium geometry. 

The Hubbard molecule allows us to look at this instability in its most simple form. 
The possibility to obtain exact results serves a guideline to improve on the current 
theory. We have demonstrated that by adding fragmentation to the single-particle 
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propagator, it is possible to go beyond the point of instability in the RPA. 

5.1 Outlook 

Although the FRPA is a promising method, there is still room for improvement, 
not only on the level of the numerical methods used, but also with respect to the 
physical content of the theory. We will list some of the possible pathways to a 
more refined Green's function method: 

• The computational cost for solving the FRPA problem can be drastically 
reduced by using Arnoldi techniques in the 2plh and 2hlp diagonalizations. 
This approach has been previously applied in Ref. [BHJ09] and it was found 
that a limited number of Arnoldi vectors guarantee correct converged results 
for total energies and ionization potentials. 

• The SCRPA (self-consistent RPA) may offer a viable alternative for the calcu- 
lation of ph and pp interactions. This scheme has the property that it solves 
the two-body problem exactly. It goes beyond the quasi-boson approximation 
by taking into account the killing condition for the excitation operators. From 
the formulation in Ref. [SS99], it can be seen that the SCRPA procedure has 
the same diagrammatical structure as the standard RPA and can also be 
translated in the diagrammatical language used to describe the FRPA. The 
occurence of single-particle density matrices in the SCRPA matrix forms an 
obstacle for the implementation. This problem can be tackled in two ways, 
either by resorting to perturbation theory of the density matrix in function of 
the SCRPA amplitudes or by using matrix relations that are system specific. 
For realistic systems, another solution has to be found. One possibility is 
to use the single-particle density matrix from a Faddeev calculation and to 
insert this into the SCRPA pp and ph equations, the solutions of which can 
then be used to build the Faddeev matrices. The matrix dimensions of this 
problem should stay the same as in the current implementation of the FRPA. 
The additional cost would be the iterative solution of the SCRPA equations, 
which should be small compared to the diagonalization in 2plh and 2hlp 
space. 

• The current approach would be more elegant if it could be translated into 
an EOM-approach along the lines of Ref. [ROW68]. However, the exist- 
ence of an excitation operator that has the current formalism as a result 
has not yet been found. In Ref. [JS11], the excitation operators as used 
for the normal RPA equations were extended to incorporate higher excita- 
tions. This was done in such a way that the excitation operator destroys 
the correlated vacuum state. One can apply this approach for the polariz- 
ation propagator to the case of the single-particle propagator by extending 
the excitation operators from the single-particle space to the 2plh and 2hlp 
space. The final matrix problem can be made Hermitian by symmetrizing 
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the anti-commutators. 
• An ultimate test for robustness of the FRPA would be to apply it to extended 
systems. For this to be possible a method must be devised to handle the 
continuous variables that are characteristic to this sort of systems. Their are 
two possible paths that were previously applied to nuclear matter[DNDW97], 
the BIN[NWR91, NWdSH93] and BAGEL[MTK88, MS93b] methods. In the 
BIN method, the energy axis is discretized and all quantities of interest are 
represented in these equidistant energy bins. In the BAGEL method, the 
Green's function is represented by a few characteristic poles, chosen in such 
a way that the lowest order moments of the exact distribution are reproduced. 
A more difficult path is to reformulate the theory where all the quantities are 
kept in their analytical form. Success is not guaranteed when pursuing this 
path. 



Appendix A 

Uniform electron gas 



In this chapter we derive the expression for the Lindhard function of the uniform 
electron gas in the TDA. The derivation goes along the lines of the derivation 
of the RPA polarization propagator[DN05], with the exception that the backward 
propagating amplitude is not present. The expression for the Lindhard function 
is used to calculate the discrete pole which arises in the limit of vanishing total 
momentum. This so-called plasmon-pole can be derived classically and is described 
correctly by the RPA, but not by the TDA. 



A.l Polarization propagator for infinite systems 

The Bethe-Salpeter equation for the RPA and TDA polarization propagator for 
continuous systems is given by 

n S M S (p,p';Q,£) = Wn (0) (p;Q,£) + n (0) (p;Q,£) 

x £(q, p|v$g|Q, P ")n SMs (p", p'; Q, E). 
p" 

(A.l) 

Here the Q represents total momentum of the ph state, while p, p' and p" are 
relative momenta. The interaction V corresponds to a local, central interaction 
without spin dependence, wich resulting in a polarization propagator which is di- 
agonal in total momentum and in spin degrees of freedom. The superscript [ph) 
indicates that the interaction is written in a basis that emphasizes the physical 
content of the interaction: it connects ph pairs 

(a^lV^ljS- 1 ) = (a5\V\fa), (A.2) 
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where the time-reversal relation is given in Eq. (2.38). With the current assump- 
tions for the interaction, we arrive at matrix elements 

<Q,p|^£|Q,p'> - ^MQ)-nP-p')). (A.3) 

Since the momentum Q is a conserved quantity, the direct term will completely 
dominate the ph interaction at small values of Q, which is the plasmon limit we 
are interested in. As a consequence, the sum over the momentum differences can 
be carried out. This leads to a new Bethe-Salpeter equation 



r{ph) 

where the definitions 



u SMs (q,e) = n(°)(Q, J B) + n( )(Q, J B)^(Q)n S M S (Q,s), (A.4) 



IIsm s (Q,E) = i]T^II SMs (p,p';Q,£) (A.5) 
p p' 

and 

n (0) (Q,£) = i£n(°>(p;Q,i5) (A.6) 
p 

have been used. Neglecting the exchange term and including a factor two for the 
spin degeneracy, the polarization propagator in the TDA or RPA becomes 

ns=o(Q ' E) = i-2V(Q)n(°)(Q, E y (A7) 

Assuming that the polarization propagator has a discrete pole at energy E p , the 
following Lehmann representation can be proposed 

Replacing the polarization propagator in Eq. (A. 7) and searching for the pole E p 
results in the equation 

l-2V{Q)m<®{Q,E p ) =0 (A.9) 

Finding a value for the plasmon pole is now reduced to finding an expression for 
the real part of the non-interacting polarization propagator. The simplest way of 
doing this, is by using the dispersion relation 

n(°) (Q ,£) = -- P n(0) ^W 

' tt J E-E'+ir) 

(A , 0) 

7T J E - E' - 17] 



A. 2. Tamm-Dancoff approximation 
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where the second and backward propagating term is retained only in the RPA and 
not in the TDA. The imaginary part of the Lindhard function for positive energy 
values can be calculated as 

m 2 E 



~4:TTh 3 Q 



when -2- < 2 and < E < 

PF 



Qpi- 



QL 

2m 



' 8tt h 3 Q 



when ^- < 2 and 

PF 



Qpi 



'8irh 3 Q 



^-(^f-f) 2 ) 



when 2 < 

PF 



(A.ll) 



A. 2 Tamm-Dancoff approximation 

Calculating the first term of Eq. (A. 10) results in the function 

1 ^(mE-Qp F + 2mE\n 



2mE 



mE 



PF 



when -2- < 2 

Pf 



In 



8ir 2 h 3 C 



-Qpf 



PF 



when 2 < 

Pf 



2mEpf 

Q 

o 2 ) 



In 



2mE-2Qp F +Q 2 



2mE-2Qp F +Q 2 
2mE-2Qp F -Q 2 



2mE+2Qp F ~Q 2 
2mE-2Qp F -Q 2 



(A.12) 

Taking the limit for vanishing momentum and writing this expression as a series 
expansion in function of the momentum, one arrives at 

^ nE) ,^m + \% + o m ). ,A,3, 



When substituted in Eq. (A. 7) this does not have the correct Q dependence to 
cancel the behavior of the Coulomb-interaction. As a consequence the plasmon 
pole will diverge and the TDA does not reproduce the classical behavior. 



A. 3 Random Phase approximation 

For the RPA derivation, both positive and negative energy parts of the 3Tl(°' are 
needed. It suffices to see that the $sTl( )((3, E) is an even function of E, so that the 
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integral in Eq. (A. 10) can be carried out. The resulting real part of the Lindhard 
function is given by a single expression, valid for all momenta 



Ml(°\Q,E) = 



Tripp 
4n 2 H 3 

x In 



-1 



Pf 
2Q 



1 



x In 



2mE + 2Qp F - Q 2 
2mE -2Qp F - Q 2 

2mE + 2Q PF + Q 2 



( mE Q 
V Qpf 2p F 

Pf 



2Q 



1 



mE 
Qpf 



2 Vl 



2mE -2Qp F - Q 2 



(A.14) 



The backward propagating term cancels the terms that are linear in the energy. 
The low-momentum limit is given by 

& wm ^ E > = ^(l%$ + °^)- < A - 15 > 

which exactly cancels the ^ of the interaction and results in the correct classical 
plasmon pole. 



Appendix B 

Pictures 



We will introduce two pictures of quantum mechanics in which we can analyze the 
Schrodinger equations. 

B.l Schrodinger picture 

Most of the time we deal with a description of quantum mechanics in which we 
assume that the state vectors are time dependent and the operators which act on 
these states are time independent. In this language the time-dependent Schrodinger 
equation is given by 

ih^ s (t)) = tf|*s(*)>, (B.l) 

where we have explicitly made clear the picture to which the states by the subscript 
S. The only part of this equation that dependends on time is the derivative in the 
left-hand side. To solve this equation, it is enough to have the state l^sC^o)) at 
some time t > after which the behavior at any time is determined by 

!*.(*)> = exp(-^H(t-t )^\* s (t ). (B.2) 

The exponential function of the Hamiltonian is the time-evolution operator in the 
Schrodinger picture and is defined by its power series expansion in function of H. 

B.2 Heisenberg picture 

In the Heisenberg picture, the states are considered to be time-independent and 
all time dependence is assumed to be included in the operators. The relation with 
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the Schrodinger states is given by 

|* H (t)) = exp Qfft) |* s (t)> (B.3) 

Combined with the Schrodinger equation (Eq. (B.l)), this shows the time-independence 
of the Heisenberg state \^h) 

ih^\* H (t)) = -Hexp (j_H?j |* s (t)> +cxp QlTt) *ft^|* s (t)jB.4) 
= 0. (B.5) 
The operators in the Heisenberg picture are time-dependent 

(%(t)\O s \*s(t)) = (^exp^HtjOsexp^—Htj^H), (B.6) 
from which their relation to the Schrodinger picture operators follow as 



H (t) = cxp l^-Htj O s cxp \-j.Ht ) . (B.7) 
Now the operators also follow an equation of motion 

ihjO H (t) = -H cxp (^HtjO s cxp (-^Hi 

+ exp(^Htjo s Hexp^-^Ht y j (B.8) 

= [0„(t),H]. (B.9) 
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The Faddeev Random Phase Approximation is a Green's function technique that makes use of Faddeev- 
equations to couple two-particle-one-hole and two-hole-one-particle excitations to the single-particle 
spectrum. Solving these equations implies an infinite partial summation of the perturbation expansion 
of the self energy. This method goes beyond the third-order Algebraic Diagrammatic Approximation by 
treating both the particle-hole and particle-particle interactions at the Random Phase Approximation- 
level. This paper presents the first results of our calculations for diatomic molecules. 
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1. Introduction 



2. Theory 



The study of electron correlation by means of first-principle 
calculations has taken a high rise thanks to modern computer 
technology. The Green's function formalism is one of these first- 
principles methods that has successfully been applied in quantum 
chemistry [1-4]. The correlations in a many-body system are de- 
scribed by an energy dependent self energy. A particular third- 
order approximation scheme to the self energy is the Algebraic 
Diagrammatic Construction (ADC) {5] developed by Schirmer and 
coworkers. This method has proven very successful in predicting 
electron properties in molecules [6]. This method is equivalent 
with resumming all two-particle-one-hole (2plh) and two-hole- 
one-particle (2hlp) interactions up to Tamm-Dancoff (TDA) level. 
It has proven very difficult to go beyond the TDA-level [7], even 
though it is known that the Random Phase Approximation (RPA) 
should be better to describe long range correlations. The Faddeev 
Random Phase Approximation (FRPA) [8] solves this problem by 
using the Faddeev technique to include RPA-phonons in the self 
energy. This method has successfully been applied to nuclei [9] 
and atoms [10]. In the first section of this work we will give a 
short overview of the working equations for the FRPA method. In 
Section 3 we will present the numerical results for a set of di- 
atomic molecules. 
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2A. Single-particle Green's function 

The behavior of an electronic many-body system is governed by 
the Hamiltonian 



(1) 



where a a [a' a ) is the annihilation (creation) operator for a state 
with quantum numbers a and where T a ,p and V a {i, y & are the ma- 
trix elements of the one-body operator and the antisymmetrized 
two-body operator respectively. For the present study, the one- 
body operator T represents the kinetic energy and the attraction 
to the nuclei and V is the Coulomb repulsion between the elec- 
trons. The evolution of this N-body system can be described by 
the single-particle propagator [11] 

G a .Mt.t') = -J«I^K(t)a;(t')]|<). (2) 

where T[. . .] represents the time ordering operator and a a (t) and 
a£(t) are the addition and removal operators in the Heisenberg 
picture. For practical calculations, the Lehmann representation of 
the Green's function 



G a ,ff(E)= 

m>F 



(<|a a |^ +1 )« +, |a> N ) 



-E 

m < F 



-IE, 



(3) 
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The Faddeev random-phase approximation is a Green's function technique that makes use of Faddeev equations 
to couple the motion of a single electron to the two-particle-one-hole and two-hole-one-particle excitations. This 
method goes beyond the frequently used third-order algebraic diagrammatic construction method: all diagrams 
involving the exchange of phonons in the particle-hole and particle- particle channel are retained, but the 
phonons are now described at the level of the random-phase approximation, which includes ground-state 
correlations, rather than at the Tamm-Dancoff approximation level, where ground-state correlations are excluded. 
Previously applied to atoms, this paper presents results for small molecules at equilibrium geometry. 
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I. INTRODUCTION 

The study of electronic systems by means of first-principles 
calculations has taken a high rise thanks to modern computer 
technology [1-5]. The Green's function formalism [6-8] is one 
of these ab initio methods that has been successfully applied in 
quantum chemistry [9-12]. The correlations in a many-body 
system are described in terms of an electron self-energy, which 
acts as an energy-dependent potential describing the motion 
of a single electron in the many-electron system. The accuracy 
of the (in principle, exact) Green's function formalism is now 
governed by the approximation chosen for the electron self- 
energy. 

A particular third-order approximation scheme to the self- 
energy can be obtained using the algebraic diagrammatic 
construction [ADC(3)] method as developed by Schirmer 
and co-workers [13]. This method has proven to be very 
successful in predicting one-electron properties in molecules 
[14] as measured, e.g., in electron momentum spectroscopy. 
Although the ADC(3) equations were derived in a purely 
algebraic manner, they can be shown to be equivalent to 
resumming all particle-hole (ph) and particle-particle (pp) 
interactions between two-particle-one-hole (2plh) and two- 
hole-one-particle (2hlp) states up to the Tamm-Dancoff 
approximation (TDA) [7] level. This is diagrammatic ally 
equivalent to considering phonons (defined here as excitations 
in the ph and pp channels) at the TDA level, and then allowing 
the exchange of these phonons in all possible ways between 
the tree propagators describing the 2plh or 2hlp states. 

The TDA allows no ground-state correlations in the 
construction of the phonons. An improvement in this respect 
is the random-phase approximation (RPA) [15]. At least for 
nuclear systems [ 1 5] , it is known that the RPA performs better 
at describing collective behavior. Calculations for the electron 
gas also show that the RPA leads to a correct prediction of the 
plasmon pole, whereas the TDA completely fails to describe 
the plasmon spectrum. It is therefore of interest to formulate 
an analogous theory to ADC(3) that resums the ph and pp 
interactions up to the RPA level, i.e., replacing the exchange of 
TDA phonons by the exchange of their RPA counterparts, with 



the ultimate purpose of arriving at a self-energy approximation 
that is valid for both finite and extended systems. 

Going beyond the TDA level has proven to be very difficult 
[16]. The Faddeev random-phase approximation (FRPA) [17] 
has solved this problem by using the Faddeev technique [18] 
to include RPA phonons in the self-energy. The FRPA method 
has been successfully applied to both nuclei [19,20] and atoms 
[21]. It is the aim of this paper to explore the application of 
this technique to simple molecular systems. 

In the second section of this paper, we give an overview of 
the working equations for the FRPA method. In Sec. Ill, we 
present the numerical results for a set of small molecules. A 
summary is provided in Sec. IV. 

II. THEORY 
A. Single-particle Green's function 

The single-particle motion in an A^-body system is de- 
scribed by the single-particle propagator [7,8] (atomic units 
are used throughout the paper) 



(1) 



where T[. . .] represents the time-ordering operator, ty^ is 
the exact ground state, and a a (t) and a'l(t) are the addition 
and removal operators in the Heisenberg representation for an 
electron in a single-particle state a. For practical calculations, 
it is more convenient to use the Lehmann representation of the 
Green's function 



E£i+n, 



E - - E£-' 

L.,„fL 



-hi 



fa,mfL 



- >n 



(2) 



'' matthias. degroote® ugent.be 



where the ^,^ ±] represent exact eigenstates of the Hamiltonian 
with energy E^ ±l . This transition to the energy domain 
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The accuracy of the Faddeev random phase approximation (FRPA) method is tested by evaluating total and 
ionization energies in the basis-set limit. A set of light atoms up to Ar is considered. Comparisons are made with 
the results of coupled-cluster singles and doubles (CCSD), with third-order algebraic diagrammatic construction 
[ADC(3)], and with the experiment. It is seen that even for two-electron systems, He and Be 2 " 1 ", the inclusion 
of RPA effects leads to satisfactory results, and therefore it does not overcorrelate the ground state. The FRPA 
becomes progressively better for larger atomic numbers, where it gives ~5 mH more correlation energy, and it 
shifts ionization potentials by 2-10 mH with respect to the similar ADC(3) method. The ionization potentials 
from FRPA tend to reduce the discrepancies with the experiment. 
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T. INTRODUCTION 

Ab initio studies of electronic systems aim at a direct 
solution of the Schrodinger equation in terms of the underlying 
Coulomb interaction, thus avoiding phenomenological input 
[1]. Due to its favorable accuracy and the gentle scaling 
of computational requirements with increasing number of 
particles the single-reference coupled-cluster (CC) method [2] 
has become the most frequently used theory in contemporary 
investigations of molecular and atomic systems. One can thus 
calculate molecules for which full configuration interaction 
(FCI) would not be feasible. Another approach with analo- 
gous characteristics is the Green's function (GF) theory (or, 
equivalently, propagator theory) [3-6]. An early scheme based 
on this approach is the outer-valence GF (OVGF) [7,8], which 
expands ionization energies up to third order in perturbation 
theory. The OVGF is very practical and computationally 
simple, and it has found several applications to studies of 
ionization spectra [8-12]. However, it becomes inaccurate 
whenever inner- and outer-valence ionization energies (IEs) 
are subject to shake-up contaminations [13-15]. In such 
cases one needs to resort at least to the third-order algebraic 
diagrammatic construction [ADC (3)] method [16,17]. The 
ADC(n) approach is an intermediate -state representation 
[18,19] of the GF made to be consistent with perturbation 
theory up to order n. Thus, it is size consistent and can 
be systematically improved by going to higher orders. For 
the one-body propagator, ADC(3) implies performing explicit 
configuration mixing between valence electrons and shake-up 
configurations, such as two-hole-one-particle (2hlp) and/or 
two-particle-one-hole (2plh) configurations [20]. These 
states are mixed together by ADC(3) theory in a Tamm- 
Dancoff approximation (TDA) fashion. The accuracy of the 
ADC(3) approach has been tested in several studies for both 
IEs [20,21] and excited states [22,23]. 
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Green's function theory has also extensive applications 
to solid-state physics, where the most successful scheme is 
the GW approximation (GWA) [24-26]. The GWA describes 
the modification of electrons through repeated interactions 
with collective particle-hole (ph) excitations of the system, 
which are described in the random phase approximation 
(RPA). The RPA is essential for extended systems because 
it screens the Coulomb interaction at large distances [5,6] 
and it guarantees finite correlation energies in metals and the 
uniform electron gas [27-30]. In contrast, the TDA plasmon 
spectrum is incorrect and even diverges at small momenta. 
The GWA, however, is not always satisfactory. As an example, 
particle-particle (pp) and hole-hole (hh) configurations, which 
would be included by ADC(3) but not in GWA, are necessary 
to explain satellite structures above and below the Fermi 
surface [31]. How to include these effects efficiently in GWA 
is still being researched [32]. Conversely, the inclusion of 
RPA in atomic and molecular studies may be advantageous for 
describing long-range (van der Waals) forces and dissociation 
processes [33-35]. Thus, a practical method that combines 
ADC(3) with RPA might become beneficial to the fields above. 

In a recent publication we have considered the ab initio 
calculation of the Ne atom using Green's function theory in the 
so-called Faddeev random phase approximation (FRPA) [36], 
which was originally proposed for studies of nuclear structure 
[37-41]. This approach completely includes 2plh and 2hlp 
states in the self-energy but expands these configurations in 
terms of couplings between valence particles and (simpler) 
ph and pp (hh) excitations that are calculated using RPA. By 
calculating these excitations in TDA one would be led back to 
the ADC(3) scheme. Stated otherwise, the FRPA can be seen as 
an extension of ADC(3) that employs RPA for pairs of electron 
and hole excitations (see also Ref. [42] for a discussion). The 
inclusion of RPA opens a possible way for treating long-range 
correlations in both finite and extended electron systems on 
a equal footing. An application of this method to extended 
systems, such as the electron gas, could be pursued following 
the approach of Ref. [43]. Before this is attempted, however, 
one wishes to better assess the quality of the method for 
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